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1.  INTRODUCTION 

1.1  Background 

For  the  past  few  years,  the  Air  Force  Rocket  Propulsion  Laboratory 
(AFRPL)  has  sponsored  a  series  of  theoretical  and  experimental  programs  on  the 
retrieval  of  plume  flow-field  properties  by  analysis  of  the  infrared  radiative 
and  absorptive  properties  of  plumes.  This  report  concludes  the  participation 
in  these  programs  by  The  Aerospace  Corporation.  The  first  phase  (Ref.  1)  was 
a  study  of  the  classic  problem  of  retrieving  radial  profiles  of  gas  tempera¬ 
ture  and  concentration  in  cylindrically  symmetric,  gaseous  plumes  from  trans¬ 
verse  emission  and  absorption  (E/A)  profiles  obtained  in  a  fixed  spectral 
bandpass.  The  E/A  profiles  are  defined  in  terms  of  the  radial  profiles  of 
pressure,  temperature,  and  concentration  (pTc  profiles)  by  integral  equations 
of  radiative  transfer.  Retrieval  of  the  pTc  profiles  from  the  E/A  profiles 
involves  a  numerical  inversion  of  these  Integral  equations. 

In  this  first  phase  study,  an  inversion  procedure  was  developed  and 
incorporated  into  the  computer  code  EMABIC.  The  inversion  algorithm  is  an 
iterative  Abel  Inversion.  The  well-known  Abel  inversion  procedure  is  valid 
for  optically  thin  sources;  for  the  general  case  of  optical  thickness,  an 
iterative  procedure  is  required.  The  code  has  been  applied  to  several  syn¬ 
thetic  and  experimental  data  and  performs  satisfactorily  as  a  diagnostic  for 
most  gas-only  plume  problems.  Some  problems  occur  when  the  temperature 
profile  has  a  deep  minimum  on  the  plume  axis  or  when  the  input  E/A  profiles 
are  particularly  noisy,  even  if  they  are  adequately  smoothed.  This  is  an 


1 .  S .  J .  Young ,  Inversion  of  Plume  Radiance  and  Absorptance  Data  for 
Temperature  and  Concentration!  AFRPL-TR-78-60,  IT  Si  Air  Force  Rocket 
Propulsion  Laboratory,  Edwarde  Air  Force  Base,  California,  29  September 
1978. 
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Inherent  feature  of  Inversion,  however,  and  is  not  restricted  to  the  method  of 


inversion.  A  similar  inversion  code  has  been  developed  at  the  Arnold  Engi¬ 
neering  Development  Center  (AEDC)  (Ref.  2).  A  random  error  propagation 
routine  was  added  to  EMAB1C  so  that  retrieval  error  could  be  estimated  auto¬ 
matically  from  E/A  measurement  error  (Ref.  3). 

The  second  phase  of  study  was  the  consideration  of  multispectral  inver¬ 
sion  and  the  effects  of  particle  loading  in  tactical  motor  plumes.  In  multi¬ 
spectral  inversion,  retrieval  is  made  on  the  basis  of  how  E/A  spectra  vary  in 
wavelength  for  a  fixed  measurement  line  of  sight.  It  was  found  that  this 
inversion  scheme  is  not  applicable  in  the  infrared  on  either  a  monochromatic 
or  wide-band  spectral  scale  near  the  exit  plane  for  small  plumes  with  mild 
temperature  gradients,  such  as  those  characteristic  of  tactical  rocket 
motors.  Even  under  ideal  circumstances,  temperature  and  concentration 
retrieval  errors  up  to  30%  were  encountered.  The  failure  of  the  method  is 
caused  by  the  lack  of  spatial  resolution  inherent  in  the  inversion  weighting 
functions.  Results  of  this  study  are  reported  in  Ref.  4.  Because  this  method 
failed  for  purely  gaseous  plumes,  it  was  never  applied  to  two-phase  plumes. 
It  was  decided  to  revert  to  the  multiposition  inversion  diagnostic  of  the 
first  study  phase  and  to  pursue  its  application  to  two-phase,  tactical  rocket 
motor  plumes. 


2.  C.  C.  Llmbaugh,  W.  T.  Bertrand,  E.  L.  Kiech,  and  T.  G.  McRae,  Noxsle  Exit 
Plane  Radiation  Diagnostic  Measurements  of  the  Improved  Transtage  Liquid 
Rocket  Injector  Program,  AEDC-TR-79-29,  ARO  Inc.,  Arnold  Engineering 
Development  Center,  Arnold  Air  Force  Station,  Tennessee,  March  1980. 

3.  S.  J.  Young,  Random  Error  Propagation  Analysis  in  the  Plume  Diagnostic 
Code  EMABIC,  AFRPL-TR-81-59 ,  U.  S.  Air  Force  Rocket  Propulsion 
Laboratory,  Edwards  Air  Force  Base,  California,  July  1981. 

4.  S.  J.  Young,  Multicolor  Inversion  Diagnostic  for  Tactical  Motor  Plumes, 
AFRPL-TR-80-30,  U.  S.  Air  Force  Rocket  Propulsion  Laboratory,  Edwards  Air 
Force  Base,  California,  May  1980. 
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The  primary  goal  of  the  first  part  of  the  third  phase  was  to  determine 
quantitative  limits  of  particle  loading  in  realistic  tactical  rocket  motor 
plumes  for  which  the  gas  properties  of  the  plume  could  be  retrieved  with 
standard  E/A  Inversion  diagnostic  without  having  to  account  for  the  radiating, 
absorbing,  and  scattering  effects  of  the  particles.  The  procedure  was  to 
assume  flow-field  properties  for  two-phase  plumes  of  interest,  generate  E/A 
profiles  with  account  of  particles  using  a  single-scattering  plume  radiation 
code  ( EAPROF ,  Ref.  5),  retrieve  the  gas  properties  from  the  E/A  profiles  with 
the  gas-only  inversion  code  EMAB1C  under  the  assumption  that  the  profiles  were 
caused  by  gas  alone,  and  compare  the  retrieved  gas  properties  with  the  assumed 
properties.  Generally,  the  degree  of  particle  loading  was  treated  as  a  param¬ 
eter.  The  work  focused  on  tactical  rocket  motors  where  the  particle  loading 
level  is  small,  that  is,  to  motors  where  the  particulate  material  is  added  to 
the  fuel  only  as  a  stabilizer  (e.g.,  AI2O3),  or  to  motors  where  the  plume 
particulate  results  from  chemical  reactions  (e.g.,  carbon)  but  not  to  motors 
in  which  the  major  fuel  is  itself  a  metal.  The  procedure  was  applied  to  a 
minimum  smoke  propellant  (MSP)  model,  an  advanced  liquid  propellant  (ALP) 
model,  and  a  reduced  smoke  propellant  (RSP)  model.  The  results  are  reported 
in  Ref.  6. 

Two  important  results  were  obtained.  The  first  is  that  the  limit  of 
particle  loading  at  which  reasonable  (<  10%  error)  retrieval  results  can  be 


5.  S.  J.  Young,  Use^s  Manual  for  the  Plume  Signature  Code  EAPROF,  AFRPL-TR- 
81-08,  U.  S.  Air  Force  Rocket  Propulsion  Laboratory,  Edwards  Air  Force 
Base,  California,  January  1981. 

6.  S.  J.  Young,  Retrieval  of  Flow-Field  Gas  Temperature  and  Concentration  of 
Low-Visibility  Propellant  Rocket  Exhaust  Plumes,  AFRPL-TR-82-13 ,  lh  ST 
Air  Force  Rocket  Propulslor.  Laboratory,  Edwards  Air  Force  Base, 
California,  March  1982. 
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obtained  Is  generally  smaller  than  the  nominal  loading  level  for  the  plume. 
The  Important  Implication  of  this  result  Is  that,  to  the  extent  that  the 
systems  studied  are  typical,  E/A  diagnostics  on  plumes  generated  by  even  low 
and  reduced  smoke-type  propellants  require  some  account  of  particle  effects. 
The  second  result  Is  that  the  maximum  loading  level  for  acceptable  gas  temper¬ 
ature  retrieval  Is  much  higher  (about  an  order  of  magnitude)  than  that  for  gas 
concentration  retrieval.  Consequently,  in  applications  where  temperature 
retrieval  is  of  primary  concern,  the  use  of  gas-only  E/A  diagnostics  may  be 
justified  even  though  the  total  retrieval  results  may  be  substantially  In 
error. 

Analysis  was  also  made  of  a  procedure  In  which  first-order  corrections 
were  made  to  the  total  gas-plus-particle  E/A  profiles  by  using  particle-only 
E/A  profiles  obtained  outside  the  gas  absorption  band.  The  corrected  profiles 
provided  better  estimates  to  the  gas-only  profiles  needed  In  the  gas-only 
inversion.  The  particle  loading  limit  for  valid  use  of  this  procedure  is  the 
value  for  which  the  total  extinction  of  radiation  by  particles  over  a  full 
diameter  of  the  plume  is  about  10%.  (Note  that  If  this  condition  is  met,  then 
the  condition  is  also  met  that  the  attenuation  by  scattering  alone  over  this 
path  is  less  than  about  10%.  The  latter  condition  is  required  by  the  single¬ 
scattering  assumption  used  in  the  work.)  For  the  two  cases  (MSP  and  RSP) 
where  the  nominal  loading  limit  roughly  satisfied  this  condition,  the  use  of 
the  procedure  resulted  in  retrieved  gas  properties  that  were  accurate  to 
within  the  convergence  criteria  set  on  the  inversion.  For  the  ALP  model,  the 
nominal  loading  level  was  well  above  this  limit,  and  the  retrieved  results 
were  poor . 
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In  the  second  part  of  the  third  phase  (Refs-  7  and  8),  the  work,  begun  in 
Ref.  6  was  expanded  In  two  significant  aspects.  First,  the  restriction  of  the 
first-order,  off-band  (FOOB)  correction  procedure  that  the  plume  be  optically 
thin  to  total  extinction  was  relaxed.  This  was  accomplished  essentially  by 
using  the  full,  two-phase,  single-scattering  radiation  model  in  the  inversion 
code  EMABIC,  as  well  as  in  the  plume  radiation  code  EAPROF.  The  second  sig¬ 
nificant  expansion  was  that  retrieval  diagnostics  for  particle  properties  were 
developed.  The  previous  work  on  two-phase  plume  diagnostic  assumed  that  all 
particle  properties  (e.g.,  species,  index  of  refraction,  spatial  and  size 
distribution)  were  known.  The  new  work  defined  procedures  for  retrieving  the 
volume  cross  sections  for  particle  absorption,  scattering,  and  extinction,  as 
well  as  the  particle  scattering  phase  function  and  radial  temperature  profile. 

As  in  the  preceding  work  of  this  phase,  the  analyses  presented  were  made 
with  synthesized  E/A/S  data,  that  is,  data  computed  from  assumed  known  plume 
properties.  No  application  to  experimental  data  was  made.  Also,  since  the 
main  emphasis  of  this  work  was  on  the  feasibility  of  retrieval  and  not  its 
practical  application,  no  treatment  of  random  or  bias  error  propagation  was 
made.  Application  of  the  inversion  procedure  was  made  to  the  MSP,  ALP,  and 
RSP  models  described  in  Ref.  5.  Results  and  discussion  of  the  application  are 
reported  in  Refs.  7  and  8. 


7.  S.  J.  Young,  Retrieval  of  Flow-Field  Temperature  and  Concentration  of 
Low-Visibility  Propellant  Rocket  Exhaust  Plumes,  AFRPL-TR-82-038,  U.S. 
Air  Force  Rocket  Propulsion  Laboratory,  Edwards  Air  Force  Base, 
California,  January  1983. 

8.  S.  J.  Young,  User's  Manual  for  the  Flow-Field  Diagnostic  Code  EMABIC, 
AFRPL-TR-82-037,  U.S.  Air  Force  Rocket  Propulsion  Laboratory,  Edwards 
Air  Force  Base,  California,  Feb. jary  1983. 
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1.2  Scope  of  Present  Study 

In  the  course  of  this  third-phase  work.  It  was  suggested  that  further 
study  be  made  on  the  retrieval  of  particle  size  distribution  and  concentration 
from  the  volume  scattering  cross  section,  volume  absorption  cross  section,  and 
scattering  phase  function  obtained  with  these  diagnostics.  Subsequent  analy¬ 
sis  has  shown  that  this  approach  is  not  likely  to  succeed.  The  major  problem 
is  that  the  retrieved  cross  sections  and  phase  function  are  obtained  at  only 
one  wavelength  (the  off-band  wavelength  mentioned  above)  and,  in  addition,  at 
a  wavelength  that  is  inappropriate  for  particle  size  retrieval.  Typically, 
the  off-band  wavelength  will  be  greater  than  about  3  pm;  but,  the  appropriate 
wavelength  for  retrieval  of  size  distributions  is  about  the  same  as  the  mean 
particle  size,  and  this  is  typically  less  than  1  pm. 

Recognizing  this  limitation,  a  new  procedure  for  obtaining  size  distribu¬ 
tions  was  formulated.  As  in  the  previous  diagnostic  schemes,  scattering  data 
are  obtained  with  a  traversing  laser  apparatus.  Actually,  data  are  obtained 
at  three  laser  wavelengths  that  center  on  and  bound  the  expected  size  distri¬ 
bution.  An  Abel  inversion  is  made  of  these  data  and  these  inversion  results 
used  in  an  iterative  inversion  scheme  (Twomey,  nonlinear)  to  retrieve  the  size 
distribution.  In  its  most  general  form,  the  radial  profile  of  the  size  dis¬ 
tribution  and  particle  loading  can  be  retrieved.  A  detailed  discussion  of  the 
method  and  an  example  of  its  application  to  the  simplified  case  where  the 
laser  apparatus  is  not  scanned  (and  thus  Abel  inversion  is  not  required)  com¬ 
prise  the  work  described  in  Section  2  of  the  present  report.  The  retrieved 
size  distribution  in  this  example  reflects  some  average  size  distribution  over 
the  plume  axial  station.  An  error  propagation  analysis  for  the  example  is 
also  presented. 
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In  addition  to  providing  a  candidate  particle  size  diagnostic  for  future 
application  to  small  rocket  motor  plumes,  the  work  of  Section  2  provides  the 
foundation  for  analysing  the  data  which  have  actually  been  taken  at  AFRPL  in 
the  Scattering-Angle/Polar Ization  experiment.  These  data  consist  of  measure¬ 
ments  of  the  polarization  of  laser  radiation  which  has  been  scattered  by  plume 
particulates.  The  analysis  is  presented  in  Section  3  and  suggests  that  the 
experimental  conditions  under  which  the  data  are  currently  taken  is  inappro¬ 
priate  and  that  no  useful  information  on  particle  size  or  distribution  can  be 
obtained. 

In  Section  4,  an  analysis  is  made  of  another  on-going  AFRPL  measurement 
diagnostic — the  single-color  transmissometer  experiment.  In  these  measure¬ 
ments,  simple  transmission  measurements  of  laser  radiation  through  a  diameter 
of  the  plume  are  oad^.  With  auxiliary  measurements  made  on  the  total  mass 
loading  of  the  plume,  these  transmission  results  can  be  used  to  infer  a  "char¬ 
acteristic"  mean  particle  size  that  is  relatively  Insensitive  to  knowledge  of 
either  the  particle  size  distribution  or  index  of  refraction.  The  results  of 
this  analyis  are  quite  encouraging  and  indicate  that  the  mean  particle  size 
may  be  retrieved  with  errors  of  <  10  to  20%  under  good  experimental  condi¬ 
tions. 
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2.  SIZE  DISTRIBUTION  RETRIEVAL  FROM 
SCATTERING  DATA 


In  this  section,  a  three-color,  laser  scattering  diagnostic  is  developed 
for  the  retrieval  of  particle  size  distributions  in  low-visibility  propellant, 
solid  rocket  motor  plumes*  An  iteration  algorithm  based  on  Twomey's  nonlinear 
iteration  method  is  employed.  An  example  application  for  AI2O3  is  given  for 
synthetic  data  at  X  =  0.308,  0.6328,  and  1.06  pm  and  seven  scattering  angles 
between  0=5  and  20  deg. 

In  Section  2.1,  the  experimental  geometry  and  required  experimental 
results  are  discussed.  The  integral  equations  relating  this  result  to  the 
desired  particle  size  distribution  is  derived  in  Section  2.2,  and  some  impor¬ 
tant  aspects  of  the  equations  are  discussed  in  Section  2.3.  The  inversion 
algorithm  that  affects  the  retrieval  is  presented  in  Section  2.4.  Hypotheti¬ 
cal  test  conditions  and  retrieval  results  are  reported  in  Sections  2.5  and 
2.6,  respectively.  A  modification  to  the  method  that  results  in  better 
retrieval  results  is  discussed  in  Section  2.7,  and  results  of  the  improved 
method  are  presented  in  Section  2.8.  Finally,  an  estimate  of  the  effects  of 
experimental  error  on  the  quality  of  inversion  is  given  in  Section  2.9. 

2.1  Required  Experimental  Data 

The  geometry  of  the  scattering  problem  is  shown  in  Fig.  1.  Light  from  a 
laser  at  wavelength  \  enters  the  plume  perpendicular  to  the  plume  axis  along 
the  optical  path  s  located  a  transverse  distance  z  from  the  plume  diameter. 
The  intensity  of  light  scattered  through  an  angle  6  into  a  sensor  a  distance  D 
away  with  entrance  aperture  A  is 
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where  I  is  the  incident  laser  power,  N(s)  is  the  particle  number  density  at 
s,  and  do/dfl  is  the  differential  scattering  cross  section  at  s.  The  integra¬ 
tion  over  s  indicates  that  all  light  along  the  chord  8  scattered  at 
angle  0  contributes  to  the  signal.  Thus,  the  field  of  view  of  the  sensor 
should  be  wide  enough  to  view  the  entire  plume  diameter.  The  validity  of  Eq. 
(1)  also  requires  that  D  »  R  where  R  is  the  plume  radius,  and  that  absorption 
and  multiple  scattering  within  the  plume  may  be  neglected. 

Define 

F(8)-M-(|~)  (2) 

o 

da  (s,0) 

g(0,s)  -  N(s)  - - -  (')) 

so  that  Eq.  (1)  may  be  written  as 

F(0)  -  /  SiD  g(0,s)  ds.  (4) 

o 

Transformation  to  cylindrical  coordinates  (Ref.  1)  with  the  assumption  that  N 
and  dag/dJJ  (and  thus  g)  are  functions  of  radial  coordinate  r  only,  and  insert¬ 
ing  the  dependence  on  z  yields  the  Abel  intergral  equation 

R 

F(6,z)  -  2  /  g(0,r)  — - - f-rTy.  (5) 

2  (r2  -  z2)172 


i 
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If  the  function  F(  8,z)  is  measured  for  all  z  In  0  <  z  <  R,  then  a  standard 
Abel  Inversion  (Ref.  1)  of  this  equation  can  be  made  to  obtain  the  radial 
function 


do  (r ,  0) 

g(  6>r)  -  N(r)  — j- -  .  (6) 

This  scattering  function  result  Is  the  primary  input  for  retrieving  the  size 
distribution. 

If  absorption  does  occur  in  the  plume,  by  either  gas  or  particle  species, 
the  function  g( 0,  r)  can  still  be  retrieved  by  making  an  auxiliary  attenuation 
scan  at  0  *  0  and  employing  an  iterative  Abel  Inversion  (Ref.  7).  Single¬ 
scattering,  however,  must  still  be  obtained. 

A  significant  simplification  occurs  for  the  retrieval  of  g(0,r)  if  there 
Is  reason  to  believe  that  the  particle  loading  size  distribution  and  differen¬ 
tial  scattering  cross  section  are  constant  in  radius  (or,  if  not,  if  one  is 
only  interested  In  retrieving  the  "average"  properties  of  the  plume).  In  this 
case,  g(  0,  r)  Is  not  (or  assumed  not)  a  function  of  r  and  Eq.  (5)  can  be 
immediately  integrated  and  the  result  solved  for  g(0)  to  obtain 


g(6) 


F(8,z) 

2  (RZ  -  z 


072  ’ 


(7) 


F(  0,z)  need  only  be  measured  at  one  value  of  z  (say  z  *  0)  to  obtain 


g(9) 


(8) 
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2.2  Inversion  Integral  Equation 

From  the  function  g(rt  6),  one  can  obtain  the  size  distribution  of  par¬ 
ticles  within  the  plume  at  each  r.  In  this  section,  the  integral  equation 
defining  g(r, 8)  in  terms  of  the  size  distribution  is  derived.  Let  n(a,r)  be 
the  particle  number  density  size  distribution  at  r,  and  a  the  particle  rad¬ 
ius.  The  total  number  density  at  r  is 

00 

N(r)  -  /  n( a, r)  da.  (9) 

0 


The  function  g(r, 8)  is 


da  (r,  8)  “  2 

g(r,  0)  -  N(r)  -  -  -  /  ua  Q(8,a,r)  n(a,r)  da  (10) 

where  Q(8,a,r)  is  the  Mie  differential  scattering  efficiency  that  obtains  for 
scattering  angle  8,  particle  radius  a,  and  the  index  of  refraction  for  the 
plume  conditions  at  r.  Since  r  is  only  a  parameter  in  Eq.  (10),  we  can  unen¬ 
cumber  the  notation  by  dropping  r  and  remembering  that  the  following  analysis 
must  be  repeated  for  each  r  of  interest. 

The  following  development  is  made  in  terms  of  the  mass  loading  size 
distribution  m(a)  rather  than  n(a)  directly.  The  relationship  between  the  two 
is 


m(a) 


n(a) 


(11) 
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where  d  is  the  bulk  material  density  of  the  particulate  species.  In  terms  of 
this  function,  we  have  from  Eq.  (10) 

00 

g(9)  “  /  m(a)  da.  (12) 

0 

In  place  of  the  variable  a,  we  use  the  dimensionless  size  parameter 


Then 


(13) 


g(9) 


3  ir 
2d  X 


Q(6,x) 

x 


m(x)  dx 


(14) 


where  m(x)  is  the  size  distribution  in  x  and  is  related  to  m(a)  by 


m(a)  -  ~ 


(15) 


Equation  (14)  is  the  final  relationship  between  the  function  g(0)  and  the  size 
distribution  function  m(x).  It  may  be  rewritten  in  the  standard  inversion 
form  (Ref.  9)  as 


g(  6)  *  /  K(6,x)  m(x)  dx  (16a) 

0 

K(0,x)  ■  2nKo(e,x>  (16b) 


9.  S.  Tworaey,  Introduction  to  the  Mathematics  of  Inversion  in  Remote  Sensing 
and  Indirect  Measurements,  Elsevier,  New  York,  1977. 
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Ko(e*x)  “  (16c) 

where  K(  0,x)  is  the  so-called  kernel  function.  This  function  is  essentially  a 
weighting  function  that  determines  the  relative  contribution  that  particles 
with  radius  a  -  Ax/2ir  make  to  the  function  g(  0  ). 

2.3  Kernel  Functions 

Equation  (16)  relates  the  measured  function  g( 9)  to  the  desired  function 
m(x)  through  an  integration  with  the  kernel  function  K(0,x).  The  ability  to 
perform  an  inversion  of  Eq.  (16),  and  the  accuracy  with  which  it  can  be  accom¬ 
plished  lies  largely  in  the  functional  form  of  the  kernel  (Ref.  9). 

Examples  of  kernels  used  in  the  present  analysis  are  shown  in  Fig.  2. 
The  function  Kq(  0,x  )  *  Q(  0,x)/x  was  evaluated  for  the  seven  angles  0  =  5.0, 
7.5,  10.0,  12.5,  15.0,  17.5,  and  20  deg  at  the  2520  values  of  x  =  0.1, 

0.2 .  252.0.  The  curves  of  Fig.  2  have  been  degraded  from  this  resolution 

to  the  constant  logarithm  increment  shown  on  the  figure.  The  calculations 
were  made  with  an  Index  of  refraction  n  *  1.75  -  0.01,  which  approximates  the 
index  for  pure  AI2O3  from  the  near  UV  to  the  near  1R  (Ref.  10).  Q(0,  x  )  was 
computed  with  standard  Mie  scattering  algorithms  (Ref.  11). 

The  0.1  resolution  on  the  x-grid  was  chosen  so  that  simple  trapezoidal 
integration  using  these  grid  points  would  yield  Integration  results  to  better 


10.  C.  B.  Ludwig  et  al.,  Standardized  Infrared  Radiation  Model  (S1RRM) , 
AFRPL-TR-81-54,  U.  S.  Air  Force  Rocket  Propulsion  Laboratory,  Edwards  Air 
Force  Base,  California,  August  1981;  American  Institute  of  Physics 
Handbook,  McGraw  Hill,  New  York,  1957. 

11.  D.  Diermend jian,  Electromagnetic  Scattering  or.  Spherical  Polydlapersions , 
Elsevier,  New  York^  1969. 
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than  ~  1%.  The  angle  grid  was  chosen  to  represent  a  possible  experimental 
array  In  which  (1)  the  spacing  between  angles  Is  large  enough  to  accommodate 
Individual  sensors  at  each  angle,  (2)  the  smallest  angle  is  far  enough  away 
from  0=0  that  no  directly  propagated  laser  radiation  is  seen,  (3)  the  scat¬ 
tering  is  largely  forward  scattering  so  that  the  largest  signals  may  be 
obtained,  and  (4)  the  total  range  of  6  is  small  enough  that  the  expected 
signal  does  not  range  over  more  than  about  an  order  of  magnitude- 

three  features  of  the  kernels  shown  in  Fig.  2  warrant  discussion  with 
respect  to  inversion.  First,  the  kernels  limit  the  range  of  x  over  which  an 
inversion  can  be  made.  The  kernel  at  0  *  5  deg  is  the  widest  and  is  sensibly 
different  from  zero  only  between  x  =  0.8  and  =  100  while  the  half  width 
extends  only  from  x  =  2.0  to  x  =  24.  Outside  this  region,  the  kernel  is 
essentially  zero  and  no  contribution  is  made  to  g(  0  )  from  m(x)  regardless  of 
what  m(x)  is  [assuming  it  is  a  realistic  m(x)].  Conversely,  no  information 
about  m(x)can  be  obtained  outside  the  region.  As  a  conservative  estimate, 
using  the  half  width,  retrieval  can  be  made  over  only  about  an  order  of  magni¬ 
tude  in  size  from  x  =  2  to  x  =  20.  If  it  is  expected  that  a  size  distribution 
is  wider  than  an  order  of  magnitude,  then  this  kernel  "window"  must  be  moved 
across  the  distribution  by  varying  the  wavelength  The  x-window  remains 
fixed,  but  the  a-window  varies  directly  with  X.  For  example, 

at  X  =  0.4  im  the  a-window  is  a  =  0.13  to  13  pm  and  at  X  *  0.8  pm  it  is  a  = 
0.25  to  25  im. 

The  second  important  feature  of  the  kernels  is  that  their  width  varies 
with  0.  As  the  angle  increases,  the  width  decreases  whereas  the  peak  height 
at  x  =  2.5  remains  reasonably  constant.  Thus,  with  changing  0,  the  kernel 
acts  as  a  variable-width  bandpass  filter  with  the  lower  cut-off  fixed.  It  is, 
in  fact,  this  feature  that  makes  an  inversion  possible  at  all.  If  the  kernels 
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did  not  change  with  6,  then  measurements  at  all  angles  would  give  the  same 
result  for  g( 6)  and  the  problem  would  be  degenerate* 

The  third  feature  worth  noting  Is  the  highly  fluctuating  structure  of  the 
kernels  even  in  the  degraded  plots  of  Fig.  2.  The  consequence  of  these  fluc¬ 
tuations  will  be  discussed  later. 

2.4  Inversion  Algorithm 

The  preceding  formulation  and  discussion  on  kernel  functions  is  relevant 
to  any  method  of  inverting  Eq.  (16)  for  the  unknown  function  m(x),  and  several 
methods  are  discussed  in  Ref.  9.  The  two  approaches  that  appear  most  applic¬ 
able  to  the  present  problem  are  constrained  linear  inversion  (the  Phillips- 
Twomey  method)  and  nonlinear  iteration.  In  this  report,  Twomey's  method  of 
nonlinear  iteration  is  employed.  This  method,  or  minor  variations  of  it,  has 
been  used  by  Twomey  in  retrieving  particle  size  distributions  from  screening 
measurements  (Refs.  9  and  12)  and  by  Hansen  (Ref.  13)  in  obtaining  atmospheric 
aerosol  distributions  from  optical  scattering  measurements. 

The  nonlinear  algorithm  as  used  here  is  shown  as  a  computational  flow¬ 
chart  in  Fig.  3.  One  begins  with  the  known  function  measurements  g^  -  g  (8^ 
on  the  angle  grid  6^  ( i  -  1,  2,  ...,  N) ,  the  known  kernel  functions  K^(x)  - 
K(0t,x),  the  maxima  (over  x)  of  the  kernel  functions  K^max,  and  a  first  guess 
solution  m(x).  A  loop  is  then  made  over  all  of  the  measurements.  First  the 
value  g^  is  calculated  that  would  result  if  the  first  guess  solution  were  the 
true  solution.  This  value  is  then  used  in  conjunction  with  the  observed  gj  to 


12.  S.  TVomey,  J.  Comp.  Physics  18,  188  (1975). 

13.  M.  Z.  Hansen,  Applied  Optics  3441  (1980). 
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compute  a  correction  parameter  £.  Note  that  if  the  guess  is  indeed  the  true 
solution,  then  £  =  0.  The  correction  parameter  is  then  used  with  the  normal¬ 
ized  kernel  K/Kmax  to  modify  the  first  guess  solution.  If  £  =  0,  no  modifica¬ 
tion  occurs.  The  loop  is  continued  over  all  of  the  angles  0^.  The  end  result 
is  one  cycle  of  the  iteration  process.  A  second  cycle  is  then  made  with  the 
final  m(x)  obtained  in  the  first  cycle  as  the  first  guess  for  the  second 
cycle.  The  process  is  continued  until  some  convergence  critereon  is  met. 

Some  of  the  modifications  of  the  method  that  have  been  used  by  Twomey  and 
Hansen,  and  in  some  of  the  preliminary  work  here  include  (1)  restrictions  on 
the  range  of  £,  (2)  variation  in  the  functional  form  of  kernel  correction,  and 
(3)  timing  of  correction.  Twomey  suggests  (Ref.  9)  that  it  might  be  advisable 
to  restrict  the  range  of  E,  in  order  to  slow  down  the  rate  of  convergence 
(although  he  gave  no  explicit  reason  why).  The  unrestricted  range  is 
-1  <£;<«>.  Preliminary  work  indicated  that  restriction  was  not  really  neces¬ 
sary  and  that  no  detrimental  effects  developed  with  use  of  the  whole  range. 
Nevertheless,  a  restriction  to  the  range  -0.99  100  was  imposed  in  the 
interest  of  prudency.  Hansen  (Ref.  13)  achieved  a  similar  kind  of  convergence 
damping  by  modifying  the  functional  form  of  kernel  correction.  Rather  than 
the  form  1  +  E;  Ki/Kiniax,  he  used  l  +  sj  £  Ki/Kimax|R  where  s  is  +1  or  -1 
depending  on  whether  £  is  +  or  -  and  R  is  an  exponent  less  than  or  equal  to 
one.  His  work  did  show  significant  changes  in  convergence  behavior  with  R.  R 
■  0.7  provided  the  stabilest  convergence.  The  third  modification,  that  is, 
the  timing  of  corrections  has  the  biggest  effect  on  the  present  work.  In 
Twomey's  original  method  of  application,  the  individual  corrections  were  not 
applied  immediately.  The  calculation  of  g|  in  the  first  step  of  the  loop  in 


Fig.  3  was  always  made  with  the  first  guess  solution,  not  the  solution 
obtained  by  modification  In  previous  loops.  The  effect  of  this  variation  of 
the  method  in  the  present  work  was  devastating.  No  convergence  at  all  could 
be  achieved.  This  result  is  In  direct  conflict  with  Twomey's  statement  (Ref. 
9)  that  there  Is  not  much  difference  between  the  two  methods. 

The  significant  advantages  of  this  nonlinear  Iteration  algorithm  are  that 
it  is  straightforward  to  implement,  stable,  transparent  as  to  how  it  works, 
and  mostly,  that  it  always  yields  a  positive  solution  so  long  as  the  initial 
guess  solution  is  positive.  Other  inversion  methods,  including  the  popular 
method  of  constrained  linear  inversion  can  give  unrealistic  negative  values 
for  m(x).  Also,  this  method  is  purported  to  be  better  able  than  others  to 
handle  retrieval  of  functions  with  large  dynamic  range.  The  one  serious  flaw 
of  the  method  will  be  pointed  out  later  along  with  a  method  for  reducing  its 
effect. 

2.5  Test  Inversion  Conditions 

In  order  to  test  the  inversion  algorithm,  a  model  test  case  was  estab¬ 
lished  using  AI2O3  as  the  particulate  species.  The  particle  size  and  distri¬ 
bution  parameters  used  were  based  largely  on  results  obtained  from  actual 
measurements  of  small  tactical  motors  (Ref.  14)  and  are  listed  in  Table  1. 

The  loading  and  size  distribution  of  partlcules  within  the  plume  were 
assumed  to  be  uniform  in  radius.  The  size  distribution  function  used  was 


14.  T.  D.  McCay  et  al.,  Laser  Mle  Scattering  Measurements  of  Partlcule  Size 
in  Solid  Rocket  Motor  Exhaust,  JANNAF  12th  Plume  Technology  Meeting, 
November  1980. 
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RETURN 

o(x) 


Fig.  3.  Twomey  Nonlinear  Inversion  Algorithm 
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Table  1 .  Pluae  and  Particle  Size  Distribution  Parameters 


Plume  Radius  R 

Typical  Particle  Size  a 

Mass  Loading  M 

Bulk  Density  P 

Particle  Loading  N 

Wavelength  1 

Index  of  Refraction  n 

Extinction  Cross  Section  a 

e 


Size  Distribution  Parameters 


S  cm 
O-5  m 

1 .0  x  10"®g/cm3 
3.7  g/cm3 
1.16  x  l05/cmJ 
0.6328  in  (HeNe) 
1.75  -  0.01 
4.4  x  10"8  cm2 


3 

*  -1 

6  m 

1.0  m 

0.5  m  (number  density  distribution) 
1.0m  (mass  loading  distribution) 
388.8  Mm7 
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M(a)  -  A  a 


(17) 


crf-3  g-BaY 

which  corresponds  to  the  large-size  peak  of  the  bimodal  distribution  used  in 
previous  work  (Ref.  6).  The  parameters  of  the  distribution  are  given  in  Table 
1.  Linear  and  logarithmic  plots  of  the  distribution  are  shown  in  Figs.  4a  and 

CO 

4b,  respectively.  (The  normalization  of  this  function  is  J  m  (a)  da  •  M;  the 
normalization  constant  is  A  ■  Y ) /  Y  MT[  (  crK)/ y]  where  T  is  the  gamma  func¬ 
tion;  and  the  relationship  between  mass  loading  and  particle  number  density 

i.  h  -  £  «  -jft  r  (^)/r(^).) 

The  mass  loading  was  taken  as  M  *  1.0  x  10  °  g/cm  »  which  corresponds  to 
a  particle  number  density  of  N  ■  1.16  x  lQ^/cm'*.  At  A  =  0.6328  pm  (HeNe 
laser),  the  extinction  cross  section  for  the  size  distribution  is  a  -  4.4  x 
10-8  cm^.  The  extinction  over  a  full  diameter  (2R  »  10  cm)  of  the  plume  is 
then  1  -  e  ^RNo  =  0.05,  which  indicates  that  the  test  case  easily  satisfies 
the  conditions  of  single-scattering  and  negligible  absorption. 

2.6  Results 

With  the  test  conditions  discussed  in  Section  2.5,  simulated  experimental 
results  for  g( 0  )  were  calculated  with  Eq.  (12).  The  result  is  shown  in  Fig. 
5  (A  =  0.6328  pm  ).  The  values  at  the  NA  “  7  angles  0  “  5.0,  7.5, 
10.0, ... ,20.0  were  then  used  as  input  to  the  nonlinear  inversion  algorithm. 
Retrieval  was  made  in  the  range  a  *  0.06  to  6.0  pm  with  MD  **  50  equal  divi¬ 
sions  of  Log(a).  A  constant  initial  guess  of  m(a)  “  0.25  x  10-8  g/cm across 
the  retrieval  range  was  used.  Convergence  was  deemed  complete  when  the  root- 
mean-square  (rms)  difference 
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W«'(J)  -  m(j),2 

l  l - - J 


j-1 


(18) 


between  two  successive  Iterations  (primed  and  unprlmed)  was  less  than  or  equal 
to  1%.  This  was  achieved  at  the  22nd  iteration.  The  retrieved  mass  distribu¬ 
tion  Is  shown  in  Figs.  6a  and  6b. 

The  results  elicit  both  optimism  and  pessimism.  On  the  optimistic  side, 
the  general  features  of  the  retrieval  are  as  expected.  In  the  window  region 
from  a  =  0.3  to  =  2.0  un,  the  retrieved  distribution  is  of  the  order  of 
magnitude  of  the  true  distribution,  whereas  well  outside  the  region  (a  <  0.1 
and  a  >  3.0  im  ),  the  initial  guess  distribution  has  not  been  severely  dis¬ 
turbed.  Also,  the  retrieval  window  corresponds  very  closely  to  the  window 
predicted  in  Section  2.3  on  the  basis  of  the  functional  form  of  the  kernel. 
The  discouraging  result  is  the  degree  of  fluctuation  of  the  retrieved  result 
in  the  window.  These  fluctuations  result  from  the  strong  Influence  of  the 
kernel  functions  on  the  solution.  Twomey  has  shown  (Refs.  9  and  12),  in  fact, 
that  the  nonlinear  inversion  algorithm  essentially  forms  a  solution  using  the 
kernel  functions  as  a  basis  set.  That  is,  the  solution  generated  is  of  the 
form 


NA 

m(x)  -  l  C.  K  (x).  (19) 

i  -  1  1  1 

On  the  other  hand,. the  kernel  functions  in  no  way  approximate  a  complete  basis 
set,  and  thus  it  is  not  possible  to  approximate  the  solution  accurately  in 
this  way. 
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6328  nm 


Distribution  Retrieval  Results  (Log  Plot) 
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Attempts  were  made  to  supress  the  fluctuation  by  implicit  (using  smaller 
values  of  MD)  and  explicit  smoothing.  In  no  case  was  significant  improvement 
achieved.  When  smaller  values  of  MD  were  used,  the  frequency  of  fluctuations 
was  reduced,  but  the  magnitude  of  fluctuations  was  not.  In  Fig.  7  an  example 
of  explicit  smoothing  of  the  form 


1  J  +  MS 

■<J>  •  r®Vr  ,  l  ,  ■<“>  <20> 


with  MS  3  3  is  shown.  Again,  the  frequency  of  fluctuations  is  reduced  but  the 
overall  agreement  between  the  retrieved  and  true  distributions  is  not  good.  A 
method  which  does  successfully  eliminate  the  strong  influence  of  the  kernels 
is  presented  in  the  next  section. 

The  convergence  stability  result  for  the  inversion  is  shown  in  Fig.  8. 
The  residual  rms  difference  between  the  reconstructed  (primed)  and  input 
(unprlmed)  profiles 


Ag 


NA  ....  ...  2 

y  ft  <L>  -MiK) 

i ,  1  gd)  J 


(21) 


is  plotted  versus  iteration.  At  convergence  (Af  <  0.01),  the  residual 
is  Ag  =2%.  This  value  is  consistent  with  Twomey's  observation  (Refs.  9  and 
12)  that  even  for  error  free  data,  the  accuracy  of  the  retrieved  solution  is 
seldom  better  than  about  IX. 
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Fig.  8.  Convergence  Residual 
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2 .7  Modified  Inversion  Method 

The  key  to  Improving  the  accuracy  of  retrieval  Is  to  use  an  Initial, 
first-guess  solution  that  Is  closer  to  the  true  solution.  According  to  Twomey 
(Refs.  9  and  12)  and  Hansen  (Ref.  13),  the  choice  of  initial  guess  does  not 
have  much  Influence  on  retrieval  results.  Some  preliminary  work  here  supports 
the  opposite  conclusion.  A  test  case  was  run  with  the  flat  size  distribution 
m(a)  =  1.0  x  10'6  g/cm^.  Synthetic  data  were  generated  and  inversion  runs 
were  made  with  constant  initial  guess  solutions  of  m(a)  =  1.0  x  10  ^ ,  0.5  x 
10-6,  and  0.1  x  10-^  g/cra^.  The  retrieved  results  are  shown  in  Fig.  9.  When 
the  Initial  guess  was  In  fact  the  true  solution,  the  true  solution  was 
retrieved.  As  the  Initial  guess  was  reduced  in  value,  the  retrieved  solution 
became  progressively  worse  with  respect  to  the  magnitude  of  fluctuations. 
These  results  support  the  conclusion  that  the  better  the  initial  guess,  the 
better  the  solution. 

In  synthetic  analyses  such  as  performed  here,  one  knows  how  to  construct 
a  "better"  first  guess  solution  because  one  knows  the  identify  of  the  true 
solution.  In  application  to  experimental  results,  one  may  also  have  some  idea 
from  auxiliary  observations  as  to  what  the  distribution  is  and  one  should 
certainly  take  full  advantage  of  the  information.  On  the  other  hand,  from  a 
purely  objective  point  of  view,  the  solution  is  completely  unknown  and  the  use 
of  some  constant  initial  guess  (as  is  done  here)  is  probably  all  that  is  jus¬ 
tified.  The  question  then  is,  how  does  one  construct  a  better  first  guess. 
The  answer  Is,  do  an  inversion  using  a  constant  first  guess  and  use  the  result 
as  the  first  guess  to  another  inversion.  After  all,  this  result  obtained  In 
Fig.  6  does  look  something  like  the  true  solution.  But,  a  direct  application 
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Uniform  Size  Distribution 


of  this  approach  cannot  work.  The  solution  shown  in  Fig.  6  is  the  result  of 
22  cycles  and  is  less  than  1%  different  from  the  result  obtained  at  the  21st 


cycle.  If  it  were  used  as  the  first  guess  solution  in  a  new  inversion,  the 
new  inversion  would  only  require  one  cycle  for  convergence  and  the  solution 
would  be  only  1%  (probably  less)  different  from  the  initial  guess.  In  effect, 
the  second  inversion  would  merely  perform  the  23rd  cycle  of  the  first  inver¬ 
sion.  What  is  needed  as  a  first  guess  for  the  second  inversion  is  something 
that  looks  like  the  solution  of  Fig.  6  but  is  not  exactly  the  same.  The 
function  used  here  is  the  smoothed  version  of  the  solution  shown  in  Fig.  7. 

When  this  smoothed  version  of  the  solution  was  used  in  a  second  inversion,  a 

marked  improvement  in  the  quality  of  the  solution  was  obtained.  But  why  stop 
here?  What  happens  if  this  solution  is  now  smoothed  and  used  as  the  initial 
guess  in  a  third  inversion?  The  answer  is  the  the  solution  is  even  better. 
In  effect,  a  super-iteration  is  performed.  And  the  process  can  be  continued 
indefinitely.  The  result  after  five  cycles  of  this  process  is  shown  in  Figs. 

10a  and  10b.  The  improvement  in  results  over  a  single  inversion  cycle  is 

substantial.  No  significant  improvement  was  achieved  beyond  five  cycles. 

2.8  Final  Results 

The  preceding  inversion  at  X  m  0.6328  ym  was  able  to  retrieve  the  center 
of  the  size  distribution  with  significant  accuracy.  No  information  was 
retrieved  in  the  wings,  however,  because  they  were  not  in  the  retrieval  win¬ 
dow.  In  order  to  effect  inversions  in  the  wings,  the  window  was  moved  by 
changing  the  wavelength.  In  order  to  retrieve  in  the  small-size  wing,  X  ■ 
0.308  un  was  used.  This  wavelength  is  achieved  with  a  XeCl  eximer  pulsed 
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Fig  10a.  Size  Distribution  Retrieval  Results  at  Fifth  Reiteration  (Linear  Plot) 


laser.  A  continuous  argon  ion  laser  operating  In  X  *  0.33  to  0.36  pa  aay  be 
more  experimentally  practical.  For  the  large-size  wing,  X  -  1.06  pa  (NdYAG 
laser)  was  choosen.  The  Index  of  refraction  at  these  two  added  wavelengths 
was  assumed  to  be  the  same  as  at  X  *  0.6328  pm.  Also,  It  was  confirmed  that 
the  conditions  of  single-scattering  and  negligible  absorption  within  the  pluae 
prevailed  at  the  new  wavelengths.  The  results  obtained  for  all  three  wave¬ 
lengths  at  the  10th  cycle  of  the  superiteration  scheme  are  shown  in  Figs.  11 
through  13. 

The  results  obtained  at  X  9  0.308  pa  are  not  as  good  as  were  obtained 
at  X  =  0.6328  pa.  First,  the  fluctuations  caused  by  the  influence  of  the 
kernel  functions  has  not  been  eliminated  as  well  as  before.  Second,  although 
the  scheme  Is  effective  in  the  small-size  region,  the  retrieved  result  Is  sig¬ 
nificantly  below  the  true  result.  The  results  obtained  at  X  •  1.06  pm  are 
also  somewhat  degraded  from  the  quality  of  the  0.6328  pm  results.  Here,  the 
scheme  is  effective  in  the  large-size  region,  but  yields  an  overpredlctlon  of 
the  distribution.  Also,  retrieval  in  the  center  portion  of  the  distribution 
is  underpredicted. 

In  Fig.  14,  the  log-plot  results  for  all  three  inversions  are  plotted 
without  the  true  distribution.  These  curves  would  represent  inversion  results 
from  experimented  data.  Without  knowledge  of  the  true  distribution,  the  only 
thing  we  know  is  that  the  X  *  0.308  pa  curve  is  probably  the  most  accurate 
for  small  sizes,  the  X  “  0.6328  pm  curve  is  probably  most  accurate  for  middle 
sizes,  and  the  X  ■  1.06  pa  curve  is  probably  'xist  accurate  for  large  sizes. 
With  just  this  information,  the  '’objective"  final  distribution  of  Figs.  15a 
and  15b  was  constructed. 


44 


ze  Distribution  Retrieval  Results  for  X  -  0.308  (in  (Log  Plot) 


Fig. 12a.  Size  Distribution  Retrieval  Results  for  X  =  0.6328  jxm  (Linear  Plot) 


.ze  Distribution  Retrieval  Results  for  \  =  1.06  fim  (Linear  Plot) 
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Now  for  the  hard  part — does  the  result  of  Fig.  15  Indicate  a  good 
retrieval  or  not?  From  the  log  plot  of  Fig.  15a,  objection  could  be  made  to 
the  underprediction  of  the  small-slze  retrieval  and  the  overprediction  of  the 
large-size  retrieval.  It  should  be  noted,  however,  that  truly  drastic  differ¬ 
ences  from  the  true  distribution  do  not  set  in  until  the  distribution  has 
dropped  almost  an  order  of  magnitude  from  its  peak  value.  On  the  other  hand, 
and  from  a  somewhat  less  objective  point  of  view,  there  are  probably  very  few 
people  familiar  with  size  distribution  retrieval  in  hostile  environments  who 
would  judge  the  quality  of  retrieval  shown  in  the  linear  plot  as  less  than 
outstanding. 

2.9  Error  Analysis 

In  this  section,  results  are  reported  for  three  studies  made  to  assess 
the  effects  of  experimental  data  error  on  the  quality  of  retrieval  using  the 
modified  Twomey  nonlinear  inversion  procedure.  The  test  case  described  above 
and  the  error-free  synthetic  data  and  retrieval  results  for  the  case  are  the 
same  as  shown  in  Figs.  5  and  15,  respectively. 

In  the  first  of  these  studies,  random  errors  were  superimposed  onto  the 
synthetic  data  of  Fig.  5  at  the  seven  angles  6  -  5.0,  7.b,  10.0,  12.5,  15.0, 
17.5,  and  20.0  deg.  A  computer  algorithm  for  generating  random  normal  combers 
was  used.  The  rms  value  of  the  superimposed  error  for  each  point  was  set  at 
5%  of  the  value  of  the  date,  point  (not  5%  of  the  peak  measured  values  as  has 
sometimes  been  done  In  the  past).  These  data  were  then  usee  iiv  the  modified 
Twomey  inversion  procedure  tc  get  the  results  of  Fig.  16.  2his  figure  is  a 
composite  of  results  for  all  three  wavelengths  (A  *  0.308,  0.6328,  and 


Fig. 16.  Composite  of  Size  Distribution  Retrieval  Results  for  5% 
Error  Superposition 


1.06  10  )•  Using  the  "objective"  assessment  that  the  0.308  pm  results  are 
probably  best  for  the  small  sizes,  the  0.6328  pm  results  are  best  for  the 
middle  sizes,  and  the  1.06  yn  results  are  best  for  the  large  sizes,  the  final 
retrieval  size  distribution  of  Fig.  17  was  determined.  Comparison  of  this 
result  with  the  error-free  result  of  Fig.  15b  shows  that  5Z  random  errors  do 
not  significantly  affect  the  quality  of  retrieval. 

In  the  second  study,  the  magnitude  of  errors  was  increased  to  10Z.  The 
composite  and  final  retrieval  distributions  are  shown  in  Figs.  18  and  19, 
respectively.  Here,  it  is  apparent  that  serious  degradation  of  the  quality  of 
retrieval  has  occurred.  While  in  the  error-free  and  5Z  random  error  cases  the 
error  in  the  peak  of  the  distribution  was  only  a  few  percent,  here  the  error 
is  nearly  a  factor  of  two.  Thus,  it  would  appear  that  successful  retrieval 
can  be  achieved  only  with  measurement  error  not  much  greater  than  5Z.  In  a 
recent  publication  (Ref.  15),  Trakhovsky  et  al.  report  this  same  maximum 
allowed  experimental  error  in  their  application  of  the  Twomey  inversion 
method. 

The  above  results  were  obtained  for  uncorrelated  errors  at  the  measure¬ 
ment  angles.  If  meet  of  the  measurement  noise  is  caused  by  fluctuations 
within  the  plume,  and  the  measurements  at  the  scattering  angles  are  made 
simultaneously  (that  is,  measurements  are  made  at  the  same  time  with  an  array 
of  detectors  rather  than  with  one  angle-scanning  detector),  it  is  likely  that 
there  will  be  a  high  degree  of  correlation  between  che  measurement  errors. 


15.  E.  Trakhovsky  et  al.,  Applied  Optics  ^1^,  3005  (1982). 
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Fig.  17a.  Final  Size  Distribution  Retrieval  Results  for  5%  Error 
Superposition  (Linear  Plot) 


Superposition  (Log  Plot) 


Error  «  10Z  (Random) 


10Z  Error  Superposition 


Superposition  (Linear  Plot) 


Superposition  (Log  Plot) 


That  Is,  If  the  measured  value  at  6  •  5  deg  Is  higher  than  the  true  value, 
than  the  measurements  at  all  the  other  angles  will  also  be  high.  In  the  third 
error  study,  the  effect  of  this  kind  of  bias  error  was  Investigated  by  per¬ 
forming  Inversions  with  sets  of  data  that  were  uniformly  10Z  higher  and  lower 
than  the  nominal  true  data.  The  results  are  shown  In  Fig.  20.  The  effect  is 
that  the  error  propagates  through  the  inversion  to  give  an  error  In  the  size 
distribution  that  is  the  same  order  of  magnitude  as  the  experimental  error. 

The  conclusion  of  these  studies  is  that  measurement  errors  up  to  5Z  can 
be  accommodated  by  the  inversion  routine  regardless  of  the  nature  of  the 
errors.  Much  beyond,  and  certaintly  before  1QZ,  random  errors  destroy  the 
usefulness  of  the  method.  However,  if  the  measurement  errors  are  highly 
correlated,  experimental  errors  up  to  the  magnitude  of  the  desired  accuracy  of 
the  retrieved  result  can  be  tolerated. 
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Fig.  20a.  Comparison  of  Retrieval  Size  Distributions  with  ±10% 
Bias  Error  (Linear  Plot) 


Error  -  102  (Bias) 


Comparison  of  Retrieval  Size  Distributions  with  ilOX 
Bias  Error  (Log  Plot) 


3.  SIZE  DISTRIBUTION  FROM 
POLARIZATION/ SCATTERING  DATA 


tn  this  section,  an  analysis  Is  made  to  determine  If  the  scattering/angle 
polarization  data  currently  being  taken  at  AFRPL  can  be  used  to  obtain  Infor¬ 
mation  on  the  particle  size  and  distribution  function  In  the  plumes  of  small 
tactical  rocket  motors  burning  low-visibility  propellants.  In  Section  3.1, 
some  relevant  definitions  are  given,  and  a  briew  review  of  data  analysis  up  to 
the  point  of  this  report  Is  made.  The  possibility  of  performing  inversion  on 
the  Individual  (perpendicular  and  parallel)  component  intensities  of  the 
scattered  radiation  Is  made  In  Section  3.2.  In  Section  3.3,  considerations 
are  made  on  Inversion  of  the  polarization  data  proper.  These  considerations 
revolve  principally  around  the  Idea  of  linearizing  the  nonlinear  equation 
relating  polarization  and  size  distribution.  The  inversion  algorithm  used  to 
test  the  "linearized”  inversion  idea  is  discussed  in  Section  3.4,  and  results 
of  test  Inversions  are  reported  .n  Section  3.5.  The  principle  result  of  the 
work  Is  that  there  Is  little  hope  that  the  polarization  data  can  yield  any 
useful  information  on  particle  size  characteristics. 

3.1  Introduction  and  Review 

The  original  goal  of  the  scatterlng-angle/polarlzation  experiment  was  the 
retrieval  of  the  particle  size  distribution  from  measurements  of  the  polariza¬ 
tion  of  laser  light  scattered  from  a  well-defined  volume  of  plume  through 
angles  from  15  to  150  deg  (15,  30.  50,  90,  130  and  150  deg).  The  polarization 
axis  of  the  incident  linearly  polarized  light  is  45  deg  to  the  scattering 
plane  and  thus  simulates  unpolarized  light  when  the  mutually  perpendicular 
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components  of  the  scattered  light  are  tsken  to  be  parallel  and  perpendicular 
to  the  scattering  plane.  The  scattered  components  are  related  to  the  unknown 
size  distribution  by 

OO 

1(6)  =  1,(0)  ■  C  /  K.  ( 0,  x)  m(x)dx  (22a) 

1  x  0 

» 

!2(0)  =  Iu(9)  -  C  /  K2(0,  x)  m(x)dx  (22b) 

where  m(x)  is  the  mass  loading  size  distribution,  x  is  the  dimensionless  size 
parameter,  C  is  a  calibration  constant,  anc  Kj  2  are  kernel  functions. 
Explicitly,  In  terms  of  the  Mle  scattering  efficiencies  for  the  perpendicular 
and  parallel  components  Qi  and  Qj, 


3v  ^ 9’  x^ 

Kl(8*x)“m  - x -  (23a) 

3  n  $2^  8’  x^ 

«2(d’ JL-ir-  (23b) 

where  d  «  particle  bulk  density  and  X  *  wavelength.  In  Section  2,  polariza¬ 
tion  was  of  no  concern,  and  only  the  total  functions  K  ■*  X2  and  Q  •  + 

Q2  were  discussed  (see  Eqs.  16). 

The  polarization  of  the  scattered  radiation  is 


1.(0)  -  I2(0) 

P<  8)  “  1^0)  +  I2(  0)  ' 


(24) 
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A  related  variable,  and  the  one  which  will  mostly  be  used  here,  Is  the  ratio 


1,(8) 

R(9)  »  •  (25) 

The  relationship  between  P  and  R  is  R  *  (1  +  P)/(l  -  P)  •  Written  explicitly, 

GO 

/  (  0,  x)  m(x)dx 

R(6)  --4 -  •  (26) 

/  K.( 0,  x)  m(x)dx 
0 

Equation  (26)  provides  the  relationship  between  the  desired  particle  size 
distribution  function  m(x)  and  the  measured  dataR( 0) . 

Apparently,  the  data  were  originally  to  be  processed  with  an  inversion 
code  developed  at  AEDC  by  Curry  et.  al  (Refs.  16  and  17).  A  constrained 
linear  Inversion  (Philllps-Twomey  method)  was  to  be  applied  to  Eq .  (24)  or 
(25)  to  retrieve  m(x).  The  problem  with  this  approach  (and  which  the  authors 
recognize)  is  that  neither  Eq.  (24)  nor  (25)  is  a  linear  integral  equation 
[Eq.  (22)  is ] ;  thus,  a  direct  application  of  the  method  is  unsuitable.  Their 
main  argument  that  purports  to  make  the  raethoa  suitable  for  linear  inversion  s 
that  the  denominator  of  Eq.  (24)  or  (25)  can  be  considered  constant  while  a 
solution  for  m(x)  in  the  numerator  is  obtained.  The  m(x)  retrieved  is  then 


16.  J.  W.  L.  Lewis  et  al..  Determination  of  the  S-ie  Distribution  Function 
for  Particle  in  a  Hypersonic  Flow  Field,  AEDC-TR-7 7-101 ,  ARO  Inc.,  Arnold 
Engineering  Development  Center,  Arnold  Air  Force  Station,  Tennessee,  July 
1978  . 

17.  B.  P  Curry  et  al.,  Development  of  Mle  Scattering  Techniques  for  In-Situ 
Particle  Diagnostics  at  AEDC,  AiDC-TR-80-3 ,  ARO  Inc.,  Arnold  Engineering 
Development  Center,  Arnold  Air  Force  Station,  Tennessee,  November  1980. 
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used  to  update  the  denominator,  and  a  new  m(x)  is  obtained.  The  iteration 
continues  until  convergence  Is  achieved.  At  first  glance,  this  approach  seems 
reasonable,  but  in  fact,  the  neglect  of  the  variation  of  the  denominator  with 
m(x)  is  a  very  serious  error.  The  correct  linearization  of  the  problem  is 
treated  later  in  Section  3.3,  and  it  will  be  shown  there  that  the  correct 
kernels  to  use  are  dramatically  different  from  the  ones  Curry  et  al.,  used. 

With  the  failure  of  the  AEDC  methods  to  make  sense  of  the  experimental 
data,  the  next  analysis  step  seems  to  have  been  an  attempt  to  see  if  the 
polarization  data  were  at  least  consistent  with  mean  particle  sizes  obtained 
from  the  one-color  transmissometer  data  and  to  resolve  ambiguities  in  the 
transmlssometer  data.  It  did  appear  to  be  able  to  this  much  (Ref.  14). 

3.2  Linear  Inversion  on  Polarization  Components 

An  Initial  idea  In  the  current  work  on  analyzing  the  AFRPL  scattering 
data  was  to  determine  if  a  linear  inversion,  such  as  described  in  Section  2 , 
could  be  performed  on  the  individual  polarization  intensity  components  1^(0) 
and  I2( 0)  since  these  quantities  are  linear  in  m(x).  Before  any  inversion  is 
attempted,  however.  It  is  useful  to  investigate  the  kernel  functions  of  the 
Integral  equations  to  see  if  an  inversion  it  e.en  feasible.  The  point  of 
comparison  for  this  investigation  Is  the  kernel  function  variations  shown  In 
Fig.  2  for  the  linear  Inversion  analysis  of  Section  2.  As  discussed  previous¬ 
ly,  the  Important  feature  of  these  kernels  that  allows  for  an  inversion  is  the 
variation  of  the  width  of  the  functions  with  0.  as  the  angle  Increases  from  5 
to  20  deg,  the  width  decreases  while  the  peak  height  a.  x  =2.5  remains  rea¬ 
sonably  constant;  the  kernels  act  as  a  variable-width  bandpass  f  Iter  with  the 
lower  cut-off  fixed.  As  0  changes,  the  proportion  of  different  regions  of  x 
that  contribute  to  the  Integral  changes  in  such  a  way  that  a  retrieval  for 
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m(x)  is  nondegenerate. 

The  kernels  for  1^(0)  at  the  three  smallest  angles  in  the  AFRPL  array 
(i.e.,  15,  30  and  50  deg)  are  shown  in  Fig.  21.  Angles  greater  than  50  deg 
have  not  been  considered  since  the  signal  levels  are  very  low  and  noisy.  The 
kernels  of  Fig.  21  have  been  scaled  to  unity  at  their  peak  value.  The  corre¬ 
sponding  kernels  for  I2( 9)  or  1(0)  are  virtually  the  same.  Clearly,  the 
angles  used  in  Fig.  21  do  not  provide  for  as  pretty  a  variation  of  width 
with  0  as  is  displayed  in  Fig.  2  or  for  as  neat  an  interpretation. 

The  only  obvious  changes  with  0  in  the  kernels  that  might  lead  to  enough 
differentiation  to  allow  an  inversion  are  the  disappearance  of  the  peak  near  x 
=6.5  when  0  in  increased  from  15  deg,  and  the  small  shifting  of  the  main  peak 
to  lower  x  as  0  increases.  Without  performing  an  actual  inversion,  it  is  dif¬ 
ficult  to  know  whether  or  not  these  changes  are  enough.* 

Furthermore,  the  overall  width  of  the  kernel  functions  In  Fig.  21  do  not 
provide  for  as  large  a  retrieval  range  as  those  in  Fig.  2.  In  Section  2,  it 
was  established  that  a  very  good  estimate  of  retrieval  range  was  the  region 
between  the  points  where  the  widest  kernel  function  achieved  half  its  peak 
height.  For  Fig.  2,  this  region  is  2  <  x  <  20  and  provides  a  range  of  about  a 
factor  of  10.  From  Fig.  21,  generously  using  the  50-deg  curve  for  the  lower 
limit  and  the  15-deg  curve  for  the  upper  limit,  the  retrieval  range  is 
l.3<x<  7.0,  a  factor  of  only  ~  5.  Even  with  a  factor  of  10,  in  the  Section  2 
work  it  was  necessary  to  use  three  wavelengths  in  order  to  x>ve  the  retrieval 
range  over  the  expected  width  of  the  size  distribution.  With  the  fixed  single 


It  is  possible  to  perform  an  eigenvalue  analysis  of  kernel  functions  to 
determine  beforehand  whether  or  not  a  successful  inversion  is  likely.  No 
such  analysis  was  performed  here.  Further  details  of  the  method  can  be 
fo-nd  in  Ref.  9. 


wavelength  of  \  =  0.5145  im  in  the  AFRPL  scattering  experiment,  retrieval 
could  be  achieved  only  in  the  radii  range  0.1<a<  0.6. 

Another  problem  with  performing  a  linear  inversion  analysis  on  1^  or  I2 
or  Ij  +  I2  is  that  the  calibration  constant  C  is  unknown  (probably)  so  that  at 
best,  a  successful  inversion  would  give  only  a  relative  profile  for  m(x).  On 
the  whole  this  approach  did  not  appear  promising,  and  it  was  decided  not  to 
procede  with  it.  (In  hindsight,  it  might  have  been  worth  doing  it,  however, 
since  the  method  that  was  pursued  failed.) 

3.3  Linearization  of  Polarization  Equation 

Analysis  of  polarization  or  ratio  data  is  beset  at  the  beginning  with  two 
significant  disadvantages.  First,  an  absolute  size  distribution  cannot  come 
from  polarization  data.  If  the  distribution  is  m(a)  =  Mf(a)  where  f(a)  is  a 
relative  distribution  and  M  is  any  constant  (say  total  mass  loading)  then  sub¬ 
stitution  into  Eq.  (26)  shows  that  R( 0)  does  not  depend  at  all  on  M  because  it 
cancels  in  the  numerator  and  denominator.  Thus,  as  for  uncalibrated  scattered 
intensity  data,  only  a  relative  distribution  can  be  retrieved.  In  the  follow¬ 
ing  analysis,  this  restriction  must  be  kept  in  mind. 

The  second  disadvantage  is,  of  course,  that  neither  P(0)  nor  R(  0)  is 
linear  in  m(x),  and  a  method  of  making  an  a^  priori  analysis  by  way  of  kernel 
functions  is  not  entirely  clear.  The  analysis  here  employs  a  first-order 
expansion  of  the  nonlinear  Eq.  (26)  in  order  to  get  an  approximate  linear 
equation.  The  goal  is  to  obtain  an  expansion  of  R(0  )  to  first  order  In  m(x) 
about  a  size  distribution  mQ(x) .  If  the  Integrations  in  both  the  numerator 
and  denominator  of  Eq.  (26)  are  approximated  by  quadrature  summations  of  the 
form 
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/  K.  ,(0,  x)  a(x)  dx  »  l  k/1,2*(9)  m  Ax  (27) 

O  ’  1 

then  the  derivative  of  R( 0)  with  respect  to  the  mass  concentration  In  the  jth 
x  Interval  Ax  Is 


»( 6) 


/  K,(0,  x)  m(x)  dx 
0 


(1 


R(9)  ) 

RTS,  X)  > 


Ax 


(28) 


where 


X  (6,  x) 

R<e- x) •  <29> 

The  first-order  expansion  of  R( 9)  about  ' 9)  is  then 


R(6)=*(0)  +  l  (m 

0  j  0  J 


«  Rq(0)  +  /  R'(9,  x)  m(x)  dx 


(30) 


where  RQ(  8)  Is  obtained  from  Eq.  (26)  using  iBq(x)  for  m(x),  3R(  0)/  j  1 0  -• 
obtained  from  Eq.  (28)  using  idq(x)  for  m(x) ,  and  K'(  0,  x)  is  the  linear 
approxtmat ton  kernel 


72 


Figure  22  shows  K'( 8,  x)  for  the  AFRPL  angle  array.  01q(x)  has  been  taken  as  a 
uniform  distribution  of  unit  mass  ii1q(x)  -  1  g/cm^. 

These  linear  approximation  kernels  are  dramatically  different  from  the 
kernels  for  the  purely  linear  case  In  that  they  are  more  widely  fluctuating 
and  assume  negative  values.  The  ability  of  these  kernels  to  be  negative  is  a 
direct  consequence  of  not  assuming  the  denominator  of  Eq.  (26)  to  be  constant, 
as  in  the  Curry  method.  The  corresponding  kernel  for  that  method  is 

K.  (  e,  x: 

K'(0,  X)  -  — - - -  .  (32) 

/  K2(0,  x)  mQ  (x)  dx 

0 

If  a  constrained  linear  inversion  were  to  be  carried  out,  it  should  be  made 
with  Eqs.  (30)  and  (31)  by  assuming  some  initial  iiiq(x),  computing  Rq(x)  from 
Eq.  (26)  and  treating  R( 0  )  -  Rq( 0  )  as  the  experimental  data.  The  new  m(x) 
obtained  by  solving  Eq .  (30)  would  then  be  a  new  initial  m^(x)  for  the  next 
iteration.  Each  step  of  the  iteration  is  a  linear  inversion,  and  Eq.  (31), 
not  Eq.  (32)  is  the  proper  inversion  kernel  to  use.  Note  that  —.like  a  purely 
linear  inversion,  the  inversion  kernel  must  be  updated  with  the  new  m^x)  at 
each  Iteration  step. 

The  ability  to  effect  an  inversion  of  P(0  )  or  R(  0  )  to  get  m(x',  by  any 
method  lies  in  the  Independence  of  the  kernels  shown  in  Fig.  22.  It  is  diffi¬ 
cult  .o  tell  just  by  looking  at  them  whether  or  not  they  are;  nor  is  it  clear 
how  to  accurately  judge  a  retrieval  range.  At  the  very  best,  retrieval  can  be 
affected  only  In  the  x  range  where  the  kernels  are  significantly  different 
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Fig.  22a.  Linear  Approximation  Scattering  Kernel  (8  -  15  deg) 


0£  =0W 


Fig.  22c.  Linear  Approximation  Scattering  Kernel  (0-50  deg) 


Linear  Approximation  Scattering  Kernel  (0  =  90  deg) 


ANG=1 50. 00 


22f.  Linear  Approximation  Scattering  Kernel  (0  =  150  deg) 


from  zero.  With  this  critereon,  it  appears  that  the  range  is  something  nar¬ 
rower  than  1  <x<  10.  This  yields  a  radii  retrieval  range  of  something  narrow¬ 
er  than  0.1  <a<l(jn  at  A  =  0.5145  in. 


3 .4  Inversion  Algorithm 

Despite  the  lack  of  a_  priori  knowledge  on  whether  or  not  the  polarization 
inversion  would  succeed,  it  was  pursued  using  a  modification  of  the  nonlinear 
method  described  in  Section  2.  To  review  this  method  for  the  purely  linear 
case,  assume  we  have  an  estimate  of  m(x),  compute  the  left-hand  side  of  the 
linear  Eq.  (8)  at  angle  0,  and  find  the  computed  value  g * ( 0)  to  be  larger  than 
the  measured  value  g( 0) .  Then,  Twomey  argues,  we  want  to  decrease  m(x)  in 
those  regions  of  x  where  K(0,x)  is  largest  [this  decreases  g'(0)  the  most  with 
the  least  change  in  m(x)]  .  He  suggests  the  correction  (see  Fig.  3) 


m(x) 


m(x) 


(1  + 


K 


K(e,  *)  j 

KmaX  (0) 


where  Kraax  ( 0)  is  the  peak  of  K( 0,  x  )  and 


K 


e) 

irsy 


1. 


(33) 


The  correction  is  applied  cumulatively  (In  a  multiplicative  sense/'  for  all  0 
in  the  angle  array.  The  rational  for  this  procedure  follows  directly  from  the 
(symbolic)  derivative 


X)  4x. 


(34) 


80 


The  corresponding  symbolic  derivative  for  the  polarization  problem  follows 
from  Rq .  (28)  and  Is 


3R(  9) 
3m(x) 


Kj  (6,  x) 

QO 

/  K_( 0,  x)  ra(x)  dx 
0 


(1 


(35) 


For  the  nonlinear  problem,  then,  following  Twomey's  reasoning.  If  we  have  an 
estimate  of  m(x),  compute  the  left-hand  side  of  Eq.  (26)  at  angle  0  and  find 
the  computed  value  R'(0  )  to  be  larger  than  the  measured  value  R(  0  ) ,  we 
should  decrease  m(x)  in  those  regions  of  x  where  K1(0,  x)  is  largest  provided 
R(0,  x)  >  R'(  0,  x)  in  the  region;  otherwise,  m(x)  should  be  Increased.  An 
algorithm  which  accomplishes  this  for  the  linear  approximation  method  and 
takes  into  consideration  the  fact  the  K'(6,  x)  may  be  negative  Is 
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where  R0  is  the  experimentally  measured  ratio  and  RQ  is  the  ratio  predicted 
using  the  guessed-at  solution  for  m(x). 

3.5  Test  Inversion  and  Results 

A  test  case  was  run  in  order  to  estimate  the  retrieval  range  for  a.  A 
constant  mass  distribution  in  0.01  <  a  <20  pm  was  assumed  and  the  ratios 
1^/12  computed  for  the  AFRPL  angle  array.  The  results  are  given  in  Table  2. 
These  values  were  then  used  as  input  data  to  the  inversion  algorithm.  The 
first-guess  solution 


m(a)  »  Ca 

was  used,  and  iteration  was  continued  until  IX  convergence  in  m(a)  was  achiev¬ 
ed.  The  results  are  shown  in  Fig.  23. 

What  is  sought  is  the  range  over  which  the  retrieval  algorithm  can  force 
the  solution  to  have  a  flat  plateau  (see  dotted  curve  of  figure).  The 
retrieved  distribution  does  display  this  feature  in  the  range  0.1  <  a 
0.6  <  but  with  large  fluctuations  superimposed.  This  range  is  consistent  with 
the  range  of  x  over  which  tVia  linear  kernels  are  nonzero  (use  Fig.  22  and  X  » 
0.5145  um  ).  The  fact  that  "he  level  of  the  plateau  is  not  the  same  as  the 
true  level  is  a  consequence  of  the  inversion  not  being  able  to  obtain  the 
absolute  magnitude  of  the  distribution. 

As  a  further  test,  an  inversion  was  made  for  the  test  conditions  used  in 
Section  2.  The  true  dlstributicr.  is  the  solid  curve  of  Fig.  24.  The  sinu- 
lated  rac+o  data  is  given  in  Table  3.  A  first-guess  solution  of  m(a)  -  1.5  x 
10  g/cnr-pm  was  used,  and  iteration  was  continued  until  IX  convergence  was 
achieved.  The  retrieved  distribution  is  the  dashed  curve  of  Fig.  24.  Keeping 
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Table  2.  Synthetic  Polarization  Data  Using 
Constant  Size  Distribution. 


8(deg) 

R 

15 

1.026 

30 

1.011 

50 

0.9536 

90 

0.7458 

130 

0.6121 

150 

0.6909 
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Size  Distribution  Retrieval  tor  Uniform  Size  Distribution 
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Fig.  24a.  Size  Distribution  Retrieval  for  Size  Distribution  of  Fig.  4  (Linear  Plot.) 


Distribution  Retrieval  for  Size  Distribution  of  Fig.  4  (Log  Plot) 


Table  3.  Synthetic  Polarization  Data  Using 
Size  Distribution  of  Pig.  4. 


8(deg)  R 


15 

1.027 

30 

1.085 

50 

0.9690 

90 

0.6043 

130 

0.5580 

150 

1.091 
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in  mind  that  we  are  looking  for  regions  where  the  retrieved  distribution  has 
the  same  shape  as  the  true  distribution  and  not  necessarily  the  same  value,  we 
find  that  the  only  region  in  which  any  reasonable  retrieval  Is  obtained  is  in 
0.3  <a<  0.5  im ,  a  pitifully  small  range. 

In  summary,  the  most  information  on  particle  size  distribution  that  can 
be  obtained  from  the  APRPL  polarization/scattering  data  is  the  shape  of  the 
distribution  in  the  size  range  0.3  <a<  0.5  ia,  and  this  assumes  the  avail¬ 
ability  of  perfectly  accurate  data.  Although  the  appropriate  calculations 
have  not  been  done  It  can  reasonably  be  predicted  on  the  basis  of  much  past 
experience,  that  if  the  experimental  errors  of  the  APRPL  date  are  propagated 
through  the  inversion,  no  Information  on  the  size  distribution  can  be 
obtained . 
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4.  TRANSMISSOMETER  DIAGNOSTIC  FOR  MEAN  RADIUS 


4 • 1  Introduction 

Ariessohn,  Self,  and  Eustts  (Ref.  18)  describe  a  two-color  transmission 
measurement  technique  for  determining  a  particular  mean  size  (the  volume- to- 
surface  mean  radius)  and  total  mass  loading  of  particles  In  a  two-phase  com¬ 
bustion  flow.  Their  method  Is  essentially  an  extension  to  absorbing  particles 
of  a  method  devised  by  Dobbins  and  Jizmagian  (Refs.  19  and  20)  for  nonabsorb¬ 
ing  particles.  The  use  of  transmission  measurements  at  two  wavelengths  In 
this  method  allows  for  the  simultaneous  determination  of  the  volume-to -surface 
mean  radius  and  the  total  mass  loading.  If  the  total  mass  loading  Is  known, 
the  method  can  be  simplified  to  a  one-color  diagnostic  for  obtaining  just  the 
volume- to- surface  mean  diameter.  The  significant  feature  of  these  one-  or 
two-color  diagnostics  Is  that  Ignorance  of  the  actual  size  distribution  of  the 
particles  or  their  complex  Index  of  refraction  (within  limits)  causes  only  a 
relatively  small  error  in  the  value  of  the  retrieved  results. 

In  this  section,  an  assessment  is  made  of  the  applicability  of  these 
methods  to  the  retrieval  of  the  volume-to-surface  mean  particle  size  and  total 
mass  loading  In  two-phase,  low-vlslbllity  propellant,  solid  tactical  rocket 
motor  plumes.  Only  the  single-color  diagnostic  was  considered  to  the  point  of 
actually  developing  a  working  diagnostic  for  analysis  of  AFRPL  data.  This 


18.  Ariessohn  et  al.,  Applied  Optics  19^  3775  (1980). 

19.  Dobbins  and  Jizmagian,  J.  Opt.  Soc.  Amer.  Ji6  1345  (1966). 

20.  Dobbins  and  Jizmagian,  J.  Opt.  Soc.  Amer.  56,  1351  (1966). 
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diagnostic  la  treated  In  Section  4.3.  An  outline  of  a  procedure  for  develop¬ 
ing  a  two-color  diagnostic  Is  present  In  Section  4.3.  Both  of  these  diagnos¬ 
tics  are  applicable  only  to  the  case  where  transmission  data  are  obtained  for 
a  single  line  of  sight  through  a  full  diameter  of  the  plume.  As  such,  the 
retrieval  particle  radius  reflects  some  average  of  the  particle  radius  over 
the  cross  sectional  area  of  the  plume  at  the  axial  station  where  the  measure¬ 
ments  are  made.  In  principle,  these  diagnostics  could  be  extended  to  trans¬ 
mission  measurements  made  in  a  lateral  scan  across  the  plume  so  that,  in 
conjunction  with  an  Abel  Inversion,  the  radial  profile  of  mean  particle  radius 
could  be  retrieved. 

4 . 2  Single-Color  Diagnostic 

In  this  section,  a  review  of  the  single-color  method  for  a  fixed  line  of 
sight  Is  given,  and  the  retrieval  error  in  the  volume- to- surface  mean  ratio 
for  AI2O3  particles  caused  by  ignorance  of  the  actual  size  distribution  and 
Index  of  refraction  Is  determined.  The  basic  one-color  transmission  method  Is 
described  In  Section  4.2.1.  Error  introduced  by  ignorance  of  the  size  distri¬ 
bution  is  considered  In  Sections  4.2.2  and  4.2.3  and  by  Ignorance  of  the  com¬ 
plex  Index  of  refraction  in  Section  4.2.4.  The  final  retrieval  and  error 
results  are  presented  In  Section  4.2.5. 

4.2.1  Basic  Retrieval  Method 

The  extinction  of  a  well-collimated  beam  of  light  traversing  a  uniform 
region  of  suspended  particles  is  governed  by  the  equation 


■  L  /  *a^  Q  (A, a)  n(a)da 

0  * 


-lnt 


(38) 


where  t  Is  the  transmittance,  L  Is  the  path  length,  a  Is  a  measure  of  particle 
size,  n(a)  is  the  particle  size  distribution  function,  Qe( X,a)  is  the  total 
extinction  efficiency  (scattering  plus  absorption),  and  X  is  the  wavelength  of 
the  radiation.  In  the  present  application,  the  particles  are  assumed  to  be 
homogeneous  and  spherical  with  a  complex  index  of  refraction  m( X) 

=  n(X)-i  k(X).  The  medium  in  which  the  particles  are  suspended  is  assumed  to 
have  an  index  mQ  =  1  -1(0).  The  size  parameter  a  is  taken  as  the  particle 
radius. 

Define  the  mean  extinction  ef f iciency  Q( X) ,  the  total  mass  loading  Cm 
and  the  volume-to-surface  mean  radius  aj2>  respectively,  by 


Q<*) 


/  sl  Q  (X,a)n(a)da 
0  _ 

"  9 

/  a  n(a)da 


(39) 


C m  *  /  (d  y  x  a^)n(a)da 


(40) 


and 


/  a^  n(a)da  -j 


*32 


/  n(a)da  a2 


'41) 


where  d  is  the  bulk  density  of  the  particle  material.  Then,  the  transmission 
equation  may  be  written  as 


-lnt 


31.  r  Q(  X) 
tn  m  a32 


(«) 
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Equation  (42)  provides  a  means  for  obtaining  a  first-order  estimate  of 
particle  size  if  Cm  is  assumed  to  be  known.  If  the  size  distribution  Is 
assumed  to  be  monodlsperse  at  particle  radius  a^  then  Q(  A)  »  Qe(A,a32)  and 
Q(A)/a^2  can  be  readily  calculated  using  standard  Mie  scattering  algorithms  as 
a  function  of  332*  An  example  plot  is  shown  in  Fig.  25  for  n  »  1.75 
and  ic  =  0.  Rather  than  a^ »  the  result  has  been  prepared  using  the  indepen¬ 
dent  variable 


P32  *  2(n-l)  x32 


(43) 


where 


2lta32 
x32  ”  A 


(44) 


Is  the  so-called  dimensionless  size  parameter.  The  reason  for  using  p^2  as 
the  Independent  variable  is  discussed  later. 

From  a  measurement  of  x  and  knowledge  of  Cm  (and  L) ,  Eq.  (42)  can  be  used 
to  obtain  an  experimental  value  of  Q/ p^2  which  in  turn  can  be  used  with  Fig. 
(25)  to  infer  two  values  of  (assuming  the  measured  value  of  Q/p32  falls 
below  the  peak  of  Fig.  (25).  The  ambiguity  of  the  result  can  be  removed  by 
several  methods.  One  method  is  simply  by  having  prior  knowledge  of  the  range 
of  particle  sizes  tn  the  plume.  The  currently  employed  AFRPL  method  is  to  use 
the  polarization  results  at  90  deg  obtained  in  the  multiangle  polarization 
experiment  to  show  tha‘  le  lower  inferred  particle  size  can  be  deleted  (Ref. 
14).  In  the  second  of  the  Dobbins  -.nd  Jlzmagian  papers  (Ref.  20),  a  means  for 
removing  the  ambiguity  by  making  transmission  measurements  at  a  second  wave¬ 
length  is  discussed. 
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As  an  example  application,  consider  the  two  measurements  made  on  M2O3  by 


McCay  et  al.  (Ref.  14).  From  their  Table  I 

t  -  0.88  C  *5.5x10  ^g/cm^ 

in 

t  »  0.87  C  *  1.0  x  10  ^g/ca^ 

m 


L  *  10  cm  X  *  0.5145  via 


From  these  data,  one  Infers  (using  d  *  3.7  g/cm^) 


q/.32  - 


11. 5  UD 

6.9  in 


-1 


-1 


and,  using  Fig.  25  and  Eqs.  (43)  and  (44) 


32 


0 .11  or  0.30 in 
0.09  or  0 .34  im 


The  higher  of  the  retrieved  radii  reproduce  the  values  reported  In  Ref.  14. 
As  a  preview  to  section  4.2.4,  It  should  be  noted  that  these  values  reproduce 
the  results  of  Ref.  14  even  though  different  values  of  the  Index  of  refraction 
were  used.  Here,  n  *  1.75-(0)i  was  employed  while  in  Ref.  14,  n  ■  1.8-0. 1(1) 
was  used. 

4.2.2  Sensitivity  to  Size  Distribution 

The  usefulness  of  the  preceedlng  approach  for  other  than  monodisperslons 
follows  from  the  observation  by  Dobbins  and  Jizmagian  (Ref.  20)  that,  for  a 
fixed  value  of  a^.  the  'unction  Q^/a^  18  very  Insensitive  to  the  size 


[ 


distribution  of  particles.  Thus,  if  a  measurement  of  t  is  made,  and  if  the 
total  mass  loading  is  known,  a  determination  of  Q(A)/a^2  can  be  made,  and  from 
this  result,  a  valid  determination  of  a follows  regardless  of  the  size 
distribution  of  the  particles. 

The  accuracy  of  retrieval  using  this  method  was  investigated  by  comput¬ 
ing  Q(A)/a^2  for  the  range  of  all  possible  unimodal  and  bimodal  rectangular 
size  distributions  possessing  a  fixed  value  of  a^  •  Although  these  distribu¬ 
tion  are  not  realistic  in  shape,  they  provide  the  degree  of  variation  required 
and  allow  easier  mathematical  manipulation  than  the  more  usually  employed 
"realistic"  distributions.  The  unimodal  distribution  is 


n(a)  =  N  f (a) 


(45) 


where  f(a)  is  the  normalized  function 


f(a) 


*'<VV 

0 


\**u 

a>a 


(46) 


u 

a^  and  au  are,  respectively,  the  lower  and  upper  bounds  of  the  distribution. 
The  moments  of  the  distribution  are 


/  .  -L. 


n+1  n+1 

a.  -  aL 


u"  *L 


The  mean  radius  and  mean  volume-to-surface  radius  are  thus 


2  2 

1  ®u  '  \  1  , 

2  i~ -~T:  "  *  2  (au  + 

u  L 


(47) 


(48) 
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V 


(49) 


3  2  2  3 

a  +  a  a.  +  a  a,  +  a^ 


4  _  4 

3  au  *L  3  ”u  "u  “L  '  “u  "L 

a32  “  J  ~3 - 5  “  T 

au  "  *L 


- j - Y 

a  +  a  a.  +  a. 

U  U  L  L 


The  use  of  F.qs.  (43)  and  (47),  along  with  the  definitions 


n  -  a/a, 


32 


\  “  V1 


32 


(50) 


nu  "  au/a32 


for  scaled  radial  variables,  and 


x 


2»a 

~r 


x32 


(51) 


for  dimensionless  size  parameters,  transforms  Eq.  (39)  Into 


Q(m,  x32) 


nlJ  , 

/  n  Qe0»,  *32n)dn  . 


(52) 


The  argument  lists  of  Q  and  Qe  have  been  changed  to  reflect  the  dependence  on 
m  ”  n-lic  and  x32 . 

The  integral  In  Eq.  (52)  was  evaluated  by  trapezoidal  quadrature  to  an 
accuracy  of  0.1X.  Variation  ^ver  all  possible  size  distributions  with  a 
fixed  a32  was  accomplished  by  varying  the  lower  bound  parameter  in  steps  of 
0.1  from  0  to  1 .  The  corresponding  upper  bound  for  each  lower  bound  was 
determined  as  the  solution  of  the  cubic  equation  [obtained  from  Eq.  (49)} 
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v3  +  (t  -  \)  +  \  §  -  \  +  \2  (y  “  "  °- 


(53) 


As  n.  varies  from  0  to  1,  n  varies  from  4/3  to  1.  Thus,  the  narrowest  dis- 

L  U 

tribution  in  the  variation  is  the  monodispersion  “  nu  “  1 ,  and  the  widest 
is  n.  ■  0,  n  *  4/3.  The  evaluation  of  the  integrals  in  Eq.  (52)  is  somewhat 

L  U 

time-consuming.  Therefore,  the  calculations  were  done  once  and  saved  on  a 
computer  permanent  file.  Beside  the  variation  over  the  11  values  of  n^, 
calculations  were  performed  for  the  41  values  of  Xj2 


x32(i)  -  (10)1/10x(i-l);  i-2,  41;  x(l)  -  0.1 


and  the  18  values  of  index  of  refraction  implied  by 


(54) 


n  -  1.70,  1.75,  1.80 


k  -  0,  0.0002,  0.002,  0.02,  0.2,  0.5. 


(55) 


A  discussion  of  these  index  values  if  given  later. 

A  typical  result  for  the  vari icion  of  Q/a32  8ize  distribution  is 
shown  in  Fig.  26.  Again,  the  parameter  p32  is  used  as  the  independent  vari¬ 
able.  The  middle  curve  of  Fig.  26  is  the  value  of  Q/p32  averaged  over  the  11 
size  distributions  for  a  fixed  n  ■  1.75,  x  *  0.  The  other  two  curves  are  the 
mean  curve  plus  or  minus  one  standard  deviation,  3 .  The  maximum  error  in 
Q/p32  occurs  at  p32  *  4  and  is  of  the  order  9Z . 

If  the  example  application  made  in  the  last  section  is  reanalyzed  with 
Fig.  26,  the  inferred  radii  (upper  values)  for  the  two  transmission  cases  are 
now 
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*  0.30  ±  .01  for  x  -  0.88 

and 

a32  *  0*37  ±  0.02  for  x  ■  0.87 

Thus,  Ignorance  of  the  size  distribution  has  introduced  an  error  of  only  3  to 
6t  in  the  retrieved  radii. 

4.2.3  Extension  to  Blmodal  Size  Distribution 

The  preceding  analysis  was  also  carried  out  for  a  blmodal  particle  size 
distribution  consisting  of  two  nonoverlapping  rectangular  size  distributions 
(Fig.  27).  The  distribution  is  specified  by  six  parameters  (the  upper  and 
lower  bound  and  relative  concentration  of  each  component) .  Two  constraints  of 
the  distribution  are  normalization 

C1  (*u  -  *L>  *  C2  <\i  '  AL>  *  1  <56> 


and  fixed  832  radius 


a32 


3  C1  (au  ~  \  )  +  C2  (Au  ~  * 

T*  ~  .  3  3.  ,  f .  3  ,  37 

C1  (au  "  *L  5  C2  (Au  ~  *L  ) 


(57) 


These  constraints  serve  to  fix  two  of  the  size  parameters.  Tn  order  to  effect 
a  variation  over  all  possible  size  distributions,  it  is  necessary  to  vary  only 
four  remaining  parameters.  Three  of  these  were  taken  as  radial  bound  ratios 


0  <  R 


0  <  p 


(58) 


(59) 
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Fig.  27.  Bimodal  Rectangular  Size  Distribution 
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I 


0  <  r 


®L 


P- 


(60) 


The  fourth  variation  is  the  relative  concentration  in  the  two  modes  and  is 
effected  by  varying  C2  from  0  to  its  maximum  value 


r  max 
C2 


1-R _ 

|  a32  (1-R)  (1-R3) 


With  C2  specified,  is  given  by 


(61) 


where 


B 


S. _  _  _1  ([B2  -  4AC)1/2-  B) 

_  max  2A 
C2 


A  -  (P3  ~  r3)  ( p-r) 

(1  -  R)  (1  -  R3) 

^2 _  1  ((  p3-r3)  (1-R)+(1-R3)  (p-r)] 

C2“aX  (l-RKl-R3) 


max 

C2 


-  P*-r* 


1-R 


(62) 


(63) 


(64) 

(65) 


In  application,  R,  p  and  r  were  varied  on  the  grid  a/Au  -  0,  0.2,  0.4,  0.6, 
0.8,  1.0  and  C2/C2,nax  on  the  grid  0,  0.2,  0.4,  0.6,  0.8,  1.0.  Given  the 
constraints  of  Eqs.  (58)  through  (60),  this  amounts  to  336  different  distribu¬ 
tions  . 

The  computation  of  Q/ P^2  for  ^36  size  distributions,  18  combinations  of 
n  and  ic  and  41  values  of  x^2  would  be  prohibitively  expensive.  Fortunately, 
the  desired  results  can  be  obtained  in  terms  of  Q  functions  computed  for  the 
unimodal  rectangular  size  distributions  and  which  are  already  computed  and 
saved.  For  the  blmodal  size  distribution,  we  have 
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Q  - 

[*(*)]  *(•)<*« 

ci(*u3-*l3> 

3 

“u 

/ 

*L 

a2  [x(a)]da 

3  a2 

3  3 

*u  "  *L 

.  C2 

3 

Au 

/ 

a2  Qe  [x(a)]da 

,7 

AuS  '  ^ 

~ 2 

where  a  is  the  Man  square 

radius 

“o  1 

3  3 

a2  -  j  IV*!!3  "  *L3)  +  C2(Au3  ~  ' 


(66) 


(67) 


In  the  first  Integral,  define 


‘32 


l  2 *  *32 


32 


-i  .  *U 

*32 


t  *L 

^  71’ 

*32 


(68) 


4  4 

4  3  *u  "*L 

*32  “  4  3  .3 

*u  *L 


and  In  the  second,  define 


a  u  2*  *32  u  _  *u 

n  “  S'  x32  X  *  nu  _  u* 


*32 


‘32 


u  *L 


*32 


(69) 


4  4 

u  3A«  -V 
*32  '  t  Ai  ;  ;  3 

U  L 
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Then,  Eq.  (66)  is  transformed  into 


Q  -  (1  -  w) 


+  w 


(\“)3  -  (03 


1,2  Qe^x32Un^n 


(70) 


where 


C,  73  37 

it  1  Cp  "  r  > 

1  +  c~ - x~ 

2  (1  -  R3) 


(71) 


From  Eq.  (52),  we  see  that  the  quantities  in  brackets  are  just  the  Q's  for  the 
unimodal  rectangular  size  distributions.  Thus 


Q(x32)  -  (1  -  W)  Qr  (x32\  n*.  n£> 


+  w  <r  <X32U*  \U«  \U> 


(72) 


where  the  subscript  r  implies  unimodal  rectangular. 

The  computation  of  Q  in  the  bimodal  size  distribution  problem  consists 
in  calculation  x32  ,  n*,  n*,  *32,  nuU  and  n^u  and  appropriately  interpolating 
on  the  saved  files.  The  two  interpolated  values  are  then  weighted  as 
indicated  in  Eq.  (72)  and  added. 

Example  results  for  Q/ p32  averaged  over  the  336  bimodal  size  distribu¬ 
tions  are  shown  in  Fig.  28.  The  middle  solid  curve  is  Q/p^;  the  upper  and 
lower  solid  curves  are  Q/p32  *  The  variation  over  the  bimodal  size  distri¬ 
butions  results  in  a  significantly  greater  a  than  did  the  averaging  over  the 
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Fig.  28a.  Mean  Scattering  Efficiency  Results  for  Bimodal  Distribution 


unimodal  distribution.  In  fact,  for  p ^  <  0.8  and  small  <,  a  Is  larger  than 
the  computed  Q/p^.  whlch  results  In  a  negative  lower  bound.  In  such  cases, 
the  lower  bound  Is  better  taken  as  the  actual  minimum  Q/p^  occurs  In  the 
variation.  This  Is  shown  In  the  figures  by  the  lower  dashed  curve.  The  upper 
dashed  curve  is  the  maximum  Q/p^  occurring  In  the  variation. 

For  fixed  n  and  k,  the  results  presented  so  far  indicated  the  degree  of 
error  to  be  expected  given  complete  Ignorance  of  the  particle  size  distribu¬ 
tion  (except  to  the  extent  that  It  is  assumed  to  be  a  unimodal  or  bimodal 
rectangular  distribution).  Note  that  the  standard  deviation  o  has  been  used 
to  measure  the  uncertainty  in  Q/p^ •  This  is  the  "reasonable”  choice.  If 
desired,  the  analysis  could  be  carried  out  for  the  "pessimistic"  worst  case 
condition  of  using  the  minimum  and  maximum  values  of  Q/p^  that  occur  in  the 
variation  with  r  (and  they  do  not  necessarily  occur  at  r  =  0  and  1). 

A. 2. 4  Sensitivity  to  Index  of  Refraction 

In  order  to  assess  the  sensitivity  of  the  retrieval  diagnostic  to 
uncertainty  in  the  complex  index  of  refraction,  it  is  necessary  to  set  some 
bounds  on  the  real  and  imaginary  parts  of  the  index.  From  the  data  compila¬ 
tion  of  the  SIRRM  report  (Ref.  10)  and  IR  Handbook  (Ref.  21),  the  range  n  - 
1.70  to  1.80  been  selected  as  representative  of  the  real  part  of  the  index  for 
AI2O3  for  all  of  the  visible  and  near  IR  spectral  regions  and  temperatures  up 
to  the  melting  point  of  AI2O3.  The  range  of  the  imaginary  part  was  selected 


21.  The  Infrared  Handbook,  Wolf  and  Zissis,  editors,  Office  of  Naval 
Research,  Department  of  the  Navy,  Washington  D.C.,  1978. 

22.  W.  L.  Konopka  et  al.,  Infrared  Optical  Properties  of  Al^Oi  Rocket 
Particles,  preprint,  Grumman  Aerospace  Corporation  Bethpage,  New  York, 
1983. 
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primarily  from  the  data  presented  by  Koaopka  et  al.  (Ref.  22)  on  their  recent 
investigations  of  shock-tube  heated  rocket  plume  AI2O3.  Although  k  ~  10-4  to 
10~7  for  pure  AI2O3,  their  findings  Indicate  that  k  may  be  as  high  as  0.01  for 
AI2O3  actually  occurring  in  plumes.  The  work  here  uses  values  up  to  k  *  0.5 
although  k  -  0.02  Is  suggested  as  a  more  realistic  upper  bound  on  the  param¬ 
eter  . 

The  variation  with  n  In  Q/p^  f°r  8  monodlsperse  distribution  is  shown 
In  Fig.  29  for  x  *■  0.  The  variation  appears  much  smaller  than  similar  varia¬ 
tions  shown  by  Ariessohn  et  al.  (Ref.  18)  because  the  scaled  variable  p^2  has 
been  used  rather  than  a^2'  The  effect  of  using  Is  place  of  X32  is  that 
most  of  the  variation  In  the  dependent  variable  Q/p  with  n  has  been  trans¬ 
ferred  into  the  Independent  variable.  The  variation  with  x  for  n  -  1.75  is 
shown  In  Fig.  30.  Significant  variations  introduced  by  Ignorance  in  both  n 
and  <  are  limited  to  the  lower  branch  of  the  Q/p  curves  below  p^  *  3.  Of  the 
two.  Ignorance  In  x  causes  the  most  variation,  particularly  for  p^  <  1.0. 

In  order  to  arrive  at  a  retrieval  error  estimate  that  reflects  Igno¬ 
rance  of  Index  of  refraction  as  well  as  ignorance  of  size  distribution,  the 
calculations  of  Section  4.2.3  were  repeated  for  all  18  combinations  of  n  * 
1.70,  1.75,  1.80  and  x  -  0,  2  x  10-4,  2  x  10'3,  2  x  10"2 ,  0.2  and  0.5.  As 
mentioned  before,  the  latter  two  values  of  x  are  probably  excessively  high  for 
AI2O3.  Even  with  carbon  contamination,  the  value  0.02  Is  a  reasonable  upper 
bound.  The  following  analysis  was  carried  out  using  both  ■  0.02  and 
0.50.  The  36  resulting  functions  Q/P32  +  a  8n^  Q/P32  ~  a  were  plotted  on  a 
common  p  axis,  and  the  upper  and  lower  envelope  of  the  curves  used  to  define 
the  final  bounds  for  Q/p^*  The  results  obtained  using  the  unlmodal  distribu¬ 
tion  are  shown  In  Fig.  31a  and  the  results  using  the  blmodal  distribution  in 
Fig.  3lb. 
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Fig.  31a.  Retrieval  Error  Bound  Caused  by  Ignorance  of  Size 

Distribution  Function  and  Index  of  Refraction  (Unimodal) 


Fig.  31b.  Retrieval  Error  Bound  Caused  by  Ignorance  of  Size 

Distribution  Function  and  Index  of  Refraction  (Bimodal) 


A  significant  feature  of  these  results  (regardless  of  size  distribution 

used)  Is  that  above  p..s  2.5,  the  error  in  Q/p,„  does  not  depend  on  ic 

32  32  max 

Below  this  value,  the  upper  bound  depends  strongly  on  tc  .At  p,,,=  0.2, 

max  3Z 

Q/p«„  approaches  a  constant  value  that  depends  on  tc  (Q/Pi«  ♦  0.6  for  ic 
3Z  max  max 

■  0.5  and  Q/p,«  +  0.025  for  tc  -  0.02).  For  the  cases  considered  here,  the 
ii  max 

ratio  between  the  two  asymptotes  is  more  than  20. 

Comparing  the  results  for  the  two  size  distributions  (Fig.  31),  two 
differences  can  be  noted.  First,  for  intermediates  p^,  the  blmodal  curve 
displays  a  somewhat  wider  error  bound  than  for  the  unimodal  case.  A  more  sig¬ 
nificant  difference  between  the  results  for  the  unimodal  and  blmodal  distribu¬ 
tions  occurs  on  the  lower  branch  between  p^  ”  0.2  and  0.8  for  the  case 

where  ic  *  0.02.  The  lower  bound  for  the  blmodal  distribution  is  higher 
max 

than  the  bound  for  the  unimodal  distribution.  The  difference  Is  caused  by  the 
use  of  Qmin  rather  than  Q  -  o  to  define  the  lower  bound  when  Q  -  a  Is 
negative. 

4.2.5  Retrieval  and  Error  Results 

Using  the  unimodal  results  of  Fig.  31,  the  retrieval  results  shown  In 
Figs.  32  through  35  were  constructed.  (Almost  exactly  the  same  results  were 
obtained  using  the  blmodal  distribution.)  Figures  32  and  34  show  p ^  ns  a 
function  of  the  experimentally  determined  value  of  Q/p^.  The  solid  curves  Itt 
these  two  figures  are  the  appropriate  ones  to  use  If  It  is  known  on  which  of 
the  two  branches  of  the  Q/p32  versus  p32  plot  one  Is  operating  (The  AFRPL 
experiment  operates  on  the  upper  branch  p32  >  3).  For  Q/p32  >1-7  (Fig.  32) 
or  0.95  (Fig.  34)  there  Is  no  distinction  between  branches,  and  the  appropri¬ 
ate  curve  to  use  is  the  dashed  curve.  The  extension  of  the  dashed  curve  below 
these  limits  Is  the  appropriate  result  to  use  if  It  Is  not  known  at  all  on 
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which  branch  one  is  operating.  Figures  33  and  35  show  the  corresponding  error 
in  P32  * 

Figure  36  shows  a  cross-plot  of  the  error  in  p^2  versus  p 32  for  both  of 

the  ic  cases.  This  result  merely  reiterates  the  previous  conclusion  that  on 
max  r 

the  lower  branch  of  the  Q/P32  curve,  uncertainty  in  the  maximum  value  of  k 
makes  a  large  difference  in  the  estimated  retrieval  error,  but  that  on  the 
upper  branch  it  is  not  so  significant.  Using  the  nominal  value  p^2  “  6  for 
the  AFRPL  experiments,  the  estimated  retrieval  errors  are  A  and  7 %  correspond¬ 
ing  to  the  assumptions,  respectively,  that  “  0.02  or  0.50. 

From  these  curves,  p ^  and  its  error  can  be  determined  from  a  experi¬ 
mentally  determined  value  of  Q/p32.  The  volume-to-surface  mean  radius  a^  is 
related  to  p ^  by  [Eqs.  (43)  and  (44)] 

P32  -  2  (n-1)  y-  a32 

and  contains  the  real  part  of  the  index  of  refraction  as  a  parameter.  Conse¬ 
quently,  the  estimated  error  in  n  has  to  be  explicitly  Included  in  order  to 
arrive  at  the  final  error  estimate  in  832*  If  is  the  error  estimate 

obtained  from  Fig.  36  and  o  is  the  estimated  error  in  n,  then  the  error  o 

n  ’a 

in  a32  Is 


For  the  present  analysis,  a  total  range  of  1.70  -  1.80  has  been  used  for  n. 

An  appropriate  o_  for  this  range  is  o  ■  (1.80  -  1.70)/  /l2  *  0.058.  Then. 

n  it  f 

for  the  AFRPL  experimental  condition  of  P32  *6,  op  “0.07  (for  ^  «  0.5), 
the  error  in  is  about  10%. 
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The  essence  of  Che  one-color  diagnostic  presented  here  has  been  coded 
Into  a  computer  program  (A32CODE)  for  analysis  of  AFRPL  data  (Ref.  23). 

4.3  Two-Color  Diagnostic 

The  basic  one-color  diagnostic  provides  for  a  retrieval  of  the  volume-to- 
surface  mean  diameter  832  if  the  mass  loading  Cm  of  the  plume  is  known.  The 
two-color  diagnostics  allow  the  simultaneous  retrieval  of  832  and  Cm.  The 
following  paragraphs  report  some  preliminary  work  that  has  done  on  the  two- 
color  diagnostic  and  outlines  how  the  work  could  be  completed  if  interest 
revives. 

To  review,  measurements  of  plume  transmittance  are  made  at  two  wave¬ 
lengths,  X^  and  Xj,  and  the  ratio 

In 

Q2  In  Tj 

determined  [see  Eq.  (42)1 .  From  plots  of  versus  2,  *32  ls  determined; 

from  plots  of  or  Q2  versus  832*  and  with  832  known,  or  Q2  is  determined; 
and  with  or  Q2  known,  Cffl  is  obtained  from  the  basic  transmission  formula 
[Eq.  (42)). 

An  Important  added  parameter  in  the  two-color  method  is  the  ratio  of  the 
two  probe  wavelengths.  The  analysis  here  has  been  carried  out  with  the  ratios 
f  -  3.16,  5.01,  and  10.  Table  4  shows  the  wavelengths  of  commonly  used 
lasers,  and  it  can  be  seen  that  many  pair  ratios  are  close  to  these  values  of 
f. 


23.  S.  J.  Young,  Description  and  Use  of  the  Single-Color  Trsnsmissometer 
Plume  Diagnostics  Code  A32CODE ,  AFRPL-TR-84 - 04§.  P.  S.  Air  Force  Rocket 
Propulsion  Laboratory,  Edwards  Air  Force  Base,  California,  August 

1984. 
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Table  4 .  Laser  Wavelengths 


X(  pm) 

Laser 

3.51 

HeXe 

3.39 

HeNe 

1.06 

NdYAG 

0.633 

NeHe 

0.515 

Argon- Ion 

0.325 

HeCd 

0.308 

HeCl 
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Results  for  the  ratio  Q1/Q2  for  the  three  values  of  f  are  shown  In  Fig. 
37.  1*  for  Xj  (X^  <  X^)*  As  for  the  •Ingle-color  net  hod,  an  average  of 
Ql/Q2  was  made  for  11  unlmodal  rectangular  size  distributions  that  cover  all 
possible  distributions  for  a  fixed  S32*  The  upper  and  lower  curves  of  these 
figures  are  +  0  at  m  *  1.70  -  0.001  and  Q^70^  -  0  at  m  *  1.80  -  0.021, 
respectively,  where  a  Is  the  standard  deviation  In  Q^/Qj  .  The  region  of 
retrieval  lies  between  the  valley  In  the  upper  curve  and  the  peak  of  the  lower 
curve.  That  Is 


0.60  <PJ2<  4.0  f  *  3.16 

0.55  <pJ2  <4.0  f  m  5.01 

0.30  <  P32  <4.0  f  -  10 

Thus ,  the  larger  the  wavelength  ratio,  the  larger  the  retrieval  range 
for  p^2<  for  f  -  10,  the  retrieval  range  covers  about  an  order  of  magnitude 
In  retrieval  radius. 

If  one  were  operating  with  f  •  10  and  had  measured  a  Q/Q  ratio  of  0.04, 
for  example,  then  the  retrieved  p ^  from  Fig.  37  would  be  p^j15  0*93  ±  0.11. 
The  error  reflects  error  caused  by  ignorance  of  size  distribution  and  index  of 
refraction  (except  that  n  Is  assumed  to  be  in  1.70  -  1.80  and  k  <0.02).  The 
value  of  p^2  and  Its  error  would  then  be  used  with  Fig.  31  (or  Figs.  32 
through  35)  to  determine  Qj  and  its  error  caused  by  ignorance  of  size  distri¬ 
bution  and  index  of  refraction.  For  example,  the  value  p^  “  0.93  +  0.11 
yields  an  upper  value  of  §1/032“  °*22  ««d  tha  value  p^  *  0.93  -  0.11  yields  a 
lower  value  of  §1/032“  0*05  to  give  a  mean  value  and  error  of  §1/032“  * 


Fig.  37b.  Error  Bounds  in 


0.08  (~  60%  error).  Use  of  this  In  the  transmission  formula  would  then  yield 
Cffl  to  an  accuracy  of  ~  60% . 

This  magnitude  of  retrieval  error  for  Cm  Is  rather  large.  One  of  the 
reasons  for  large  error  is  that  the  retrieval  for  p^2  from  Q ^ /q2  occurs  In  the 
region  of  p^2  that  Is  below  the  peak  of  the  Q/Pj2  versus  p^2  curve  (i.e.,  on 
the  lower  branch),  and  It  Is  this  region  where  the  error  in  Q/P32  largest. 
This  problem  can  be  overcome  by  scaling  the  retrieved  value  of  p^2  and  its 
error  to  Xj.  In  this  case,  since  f  *  Xj/X^  ”  10,  the  scale  value  is  p^2  *  9 .3 
±  1.1.  Then,  again  using  Fig.  31  but  using  9.3  +  1.1  and  9.3  -  1.1  to  deter¬ 
mine  the  lower  and  upper  bounds  of  0^32“  *26  +  0.06  ( ~  23%  error). 

Actually,  the  retrieval  error  in  Cm  can  probably  be  reduced  even  below 
this  value  by  a  more  careful  analysis.  Before  abandoning  this  problem,  the 
following  procedure  was  devised  which,  if  Implemented,  would  yield  a  more 
precise  diagnostic.  The  main  idea  of  the  new  procedure  is  that  the  retrieval 


of  p32  and  Cm  would  be  carried  out  for  a  specified  n,  k,  and  rectangular  sire 
distribution.  In  effect,  plots  of  p^2  and  (and  or  Q2)  versus  (with  f 
fixed)  would  be  generated.  These  plots  would  then  be  regenerated  for  all 
n,  k,  and  size  distributions.  From  this  manifold  of  plots,  the  mean  value  of 


p^2  (and  its  standard  deviation)  and  (and  its  standard  deviation)  would  be 
calculated.  The  final  result  would  be  p^2  +  a  and  pj2  -  0  on  a  single  plot 


versus  +  Oq  and  -  Op  on  a  single  plot  versus  and 

Q2  +  Op  and  Q2  ~  °q  on  a  single  plot  versus  In  application,  one 

would  measure  Q1/Q2  and  go  directly  to  the  three  plots  to  determine  p3«,  Q1/Q2 


and  their  associated  errors.  From  p^2,  one  calculates  832  as  in  the  one-color 


method.  From  and  Q2,  one  determines  two  (hopefully  consistent)  values  of 


Cn  via  the  measured  Tj,  t2,  and  transmission  formula,  Eq.  (42). 
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